{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><center> Monte Carlo Approximation</center></h1>\n",
    "# 1. Basic\n",
    "So far, we discussed various deterministic algorithms for posterior inference. The trouble with these methods is that they can be rather complicated to derive and somewhat limited in their domain of applicability.\n",
    "\n",
    "In this chapter, we discuss an alternative class of algorithms based on the idea of Monte Carlo approximation. The idea is very simple: generate some samples from the posterior, $\\theta^s \\sim p(\\theta|D)$, and then use these to compute any quantity of interest, such as posterior predictive $p(x_{new}|D)$. All of these quantities $\\mathbb{E}(f)=\\int f(\\theta)p(\\theta)d\\theta$ can be approximated by\n",
    "$$\n",
    "\\hat{\\mathit{f}} \\approx \\frac{1}{S}\\sum_{s=1}^S f(\\theta^s)\n",
    "$$\n",
    "for some suitable function $f$. As long as the samples $\\theta^s \\sim p(\\theta)$ , then $\\mathbb{E}\\left[\\hat{\\mathit{f}}\\right]=\\mathbb{E}\\left[f\\right]$ and so the estimator $\\hat{\\mathit{f}}$ has the correct mean. The variance of the estimator is given by\n",
    "$$\n",
    "var\\left[\\hat{f}\\right]=\\frac{1}{L}\\mathbb{E}\\left[(f-\\mathbb{E}\\left[f\\right])^2\\right]\n",
    "$$\n",
    "The main issue is:how do we efficiently generate samples from a probability distribution, particularly in high dimensions. We will firstly introduce the non-iterative method for generating independent samples.\n",
    "\n",
    "# 2. Sampling from standard distributions\n",
    "## 2.1  Using the cdf\n",
    "The simplest method for sampling from a univariate distribution is based on the **invese probability transform**.Let $F$ be a cdf of some distribution we want to sample from, and let $F^{-1}$ be its inverse. Then we have the following result: If $U \\sim U(0,1)$ is a uniform random variable, then $F^{-1}(U) \\sim F$.\n",
    "\n",
    "We can take the Gaussian distribution as a example.\n",
    "$$\n",
    "p(x)=\\frac{1}{\\sqrt{2\\pi}}exp\\left(\\frac{-x^2}{2}\\right)\n",
    "$$\n",
    "The cdf of the standard Gaussian distribution is:\n",
    "\\begin{align}\n",
    "\\Phi(x)&=\\frac{1}{\\sqrt{2\\pi}}\\int_{-\\infty}^x exp\\left(\\frac{-t^2}{2}\\right)dt  \\\\\n",
    "       &=\\frac{1}{2}\\left[1+erf(\\frac{x}{\\sqrt{2}})\\right]\n",
    "\\end{align}\n",
    "where the error function is given by:\n",
    "$$\n",
    "erf(x)=\\frac{2}{\\sqrt{\\pi}}\\int_0^x exp\\left(-t^2\\right)dt\n",
    "$$\n",
    "We can get the inverse function of cdf as \n",
    "$$\n",
    "\\Phi^{-1}(p)=\\sqrt{2}\\,erfInv (2p-1)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The time of sampling 200 samples using inverse cdf algorithm is 0.000104\n",
      "The time of sampling 200 samples using np.random.randn is 0.000046\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAFpCAYAAABnHGgVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4lFX6//H3SQIJEEISEnpooUsnKL2jIgLqShMQXevuqth2F0XX8tt1rfu1rOuuDVRArAgISBOUDqGXgIQeQg2hpZB2fn+cmckkJDCBzDxT7td15RrmmcnMyYR85sx9znOO0lojhBDCvwRZ3QAhhBDlT8JdCCH8kIS7EEL4IQl3IYTwQxLuQgjhhyTchRDCD0m4CyGEH5JwF0IIPyThLoQQfkjCXQgh/FCIVU8cExOjGzZsaNXTCyGET9qwYcMprXXsle5nWbg3bNiQxMREq55eCCF8klLqoCv3k7KMEEL4IQl3IYTwQxLuQgjhhyTchRDCD0m4CyGEH5JwF0IIPyThLoQQfkjCXQgh/JCEuxBC+CEJdyGE8EMS7kII4YcsW1tGCCvs2weHDkGdOtC0KShldYuEcA/puYuAsGkTdOkC8fHQty80bw49esCqVVa3TAj3kHAXfm/6dOjaFVJS4M03YdEiePttOHAAeveGb76xuoVClD8pywi/9vPPcPfdppf+zTcQa1sFe8AAuOceuPVWGDUKKlaEYcMsbaoQ5Up67sJvHTgAd94JLVrA7NmFwW5XrRr89BN06AD33guHD1vSTCHcQsJd+CWt4U9/gpwcE+wRESXfr0oV+PJLyM2FcePM9wnhDyTchV+aORPmzYP/90wmjfP3wIULpd63aVN46y345Repvwv/4VK4K6VuVkrtVkolK6UmXuZ+dyqltFIqofyaKETZ5ObCU49epG34Xh59PhKaNYPoaPjd72Dv3hK/5777oE0bmDgRLl70cIOFcIMrhrtSKhh4HxgEtAJGK6ValXC/qsBjwNrybqQQZfHFg8s5kBrKP4OfJ+S5ifD55/DYY7BwIbRtC/PnX/I9wcFmJs3+/fC//1nQaCHKmdJXKDIqpboCL2qtb7JdfwZAa/3PYvd7G1gMPA08rbW+7O7XCQkJWjbIFuUt96MptHiwJ9ER+azbH4uKjiq88cgRGDIEtm41I6kDBlzy/b16wcGDpoMfInPJhBdSSm3QWl+xOuJKWaYu4DyPIMV2zPnJOgBxWusfy9RKIcrTmjV884ef2Uc8z3/aqGiwA9StC8uWQcuWMHx4iSWap582Z7BK7V34OlfCvaQTtB3dfaVUEPB/wFNXfCClHlRKJSqlEk+ePOl6K4W4ksxMGDeO/4Q8StP4fG69vULJ94uIgDlzzLSYe++FgoIiN996q5k6+cYbMnNG+DZXwj0FiHO6Xg9IdbpeFWgNLFNKHQC6ALNLGlTVWn+otU7QWifEFp90LMS1ePFFtiZXYuXFzjz8x2CCLvc/u2FDc4rq8uWXFNiDgmDCBLNcwfr1bm2xEG7lSrivB5oqpRoppSoCo4DZ9hu11me11jFa64Za64bAGmDolWruQpSbffvg7bf5oPk7hIWZM0+vaPx46NMH/vY3OHu2yE2jR0OlSvDJJ+5orBCeccVw11rnAY8AC4Ak4Gut9Q6l1MtKqaHubqAQVzRpElnB4UxP7c2IEWbW4xUpZSa3p6XBa68VualaNVOS//JLyMhwT5OFcDeX5rlrredprZtpreO11v+wHfub1np2CfftI7124TG7dsFXX/Hjze9x7nwQ48eX4Xs7doSRI+G99+D06SI33XcfnD8P335bvs0VwlPkDFXh2157DcLCmJp9J3XqmFUey+TZZ83Zq+++W+Rwz57QuLHpvQvhiyTche86ehSmTiVt7ATmLwll9GhzMlKZtGkDQ4fC++9DdrbjsFIwYgQsXgynTpVvs4XwBAl34bv+9z/Iz+e7+k+Qmwtjxlzl40yYYBL8q6+KHB45EvLz4fvvr72pQnjaFc9QdRc5Q1Vck9xcqF8fOnZkUMFc9uyBPXuucts8rU0PvlKlIvMftTZz3uPiTA9eCG9QnmeoCuF95s+HY8c4O+4RliyB22+/hv1QlYIHHoDERNixo8jhkSNh6VKQc+6Er5FwF75p6lSIjWVe7kByc+G2267x8ewF+y++KHL4ttvMSazz5l3j4wvhYRLuwvecOWN24Bg1iplzQqhZ0+yRek1q1IBBg8ybRn6+43CHDmZJmjlzrvHxhfAwCXfhe779Fi5eJGfU3cyfbya7XHa5AVeNG2dWjly2zHFIKbOQ5IIFss678C0S7sL3TJ0KzZqxIqsTFy7A4MHl9LhDhpjTUz///JLDFy4UyXwhvJ6Eu/AtBw+a/fDGjWP+T4oKFaBfv3J67EqVzLoD331XZN2Bfv2gcmX4URa0Fj5Ewl34lq+/NpdjxjB/vjmTtGrVcnz8u+4ywf7TT45DYWHQt68pzQjhKyTchW+ZNQs6dOBwSCN27DBjoOWqZ0+z8tisWUUO33STmUe/f385P58QbiLhLnzHiROwahUMG+boWJd7uIeEmCL+3LmQl+c4fOON5nLRonJ+PiHcRMJd+I65c81po0OH8vPPULs2tLpkq/ZyMGyYWSVyxQrHoWbNzAmxUpoRvkLCXfiOWbMgLg7drj1Ll5o6+FWflXo5N90EoaHwww+OQ0qZ3vuSJUU69EJ4LQl34RuysmDhQhg6lN2/KY4fNxspuUV4OAwYYN5MnNZeGjjQbNq0caObnleIciThLnzD4sUm4IcNY+lSc6hvXzc+37BhcOAAbNvmONSrl7n89Vc3Pq8Q5UTCXfiG2bMhIgJ692bZMrMkQHy8G59vyBBz6bTuQK1a0Ly5mWYvhLeTcBfeT2sz73zgQHSFiixb5sZ6u12tWmZhmWLTY3r1guXLiyw/I4RXknAX3m/3bkhJgYEDSUoyMyLdVm93NnCgmXp54YLjUO/epu7uVK0RwitJuAvvZ+89DxzomXq73Y03mk1BnOow9rq7lGaEt5NwF95v0SKzW3XjxixbZnZGatTIA8/bvbtZe8CpNGN/bgl34e0k3IV3y801yzEOHEhBAZ6pt9uFhZk6zMKFRQ737m1mzFi0Q6UQLpFwF95t7Vo4fx4GDmTnTrOPtUfq7XYDB0JSkqn52/TqBWlpsHOnB9shRBlJuAvvtmiR2YmjXz/P1tvtBg4sbIdN797mUkozwptJuAvvtmgRJCRAVBQrV0K9etCwoQefv00bqFmzSLg3amTm2cvJTMKbSbgL73X+PKxb5+g9r15dDnullpVS0L+/KfbbiuxKmd77L79I3V14Lwl34b1WrTJnC/XpQ2oqHDpkQbiDKbIfPQrJyY5DvXvDsWNFDgnhVSTchff69VezvnrXrqxZYw5ZEu4lFNm7dTOX9nYJ4W0k3IX3+vVX6NQJqlRh9WqoWNGsCOBxzZtDjRpFiuwtW5rFIyXchbeScBfeKSvL1Nttp4SuXg0dO5pl1j1OKdMOp557cDBcf72ZqSmEN5JwF95p7VrIyYFevcjJgcREi0oydr16maL/wYOOQzfcAFu2mPchIbyNhLvwTr/+anrMPXqwZQtcvGhxuJdQd+/SxezKJJt3CG8k4S6806+/Qrt2EBnJ6tXmkKXh3ro1REUVCfcbbjCXUncX3kjCXXifnBwzDdKp3l6vnvmyTFAQ9OxZZFC1Zk1zQpWEu/BGEu7C+2zYYArZTuHepYvFbQLTnuRkM8HdpksXGVQV3knCXXgfex2me3eOHjVjmJaWZOzsk9vt7cOUZg4fhtRUi9okRCkk3IX3Wb3a1Dtq1bL25KXiOnY0k+1XrXIcsn+ikN678DYS7sL7rFnjSPM1a6BCBZOrlgsNNYuYOYV7hw4m76XuLryNhLvwLikp5svWJV6/Htq3t+jkpZJ07WrGBC5eBEy72reXcBfeR8JdeBeneY8FBSZHExKsbVIR3bqZYN+0yXGoSxdzklVenoXtEqIYCXfhXdasMdvbtWtHcjKcO+dl4W4v/juVZm64ATIzYft2i9okRAkk3IV3Wb3aLBZWsSKJieaQV4V77dpmsNdpxkznzubS3l4hvIGEu/AeFy+ac/lt9fbERNOJb9XK4nYV162b6bnbduqIj4eICFNCEsJbSLgL77F5c5FFZDZsMIOVISEWt6u4bt1w7B6COXm1Y0cJd+FdJNyF97BPOenShfx804n3qpKMnb3u7lSa6dQJtm6F3FyL2iREMRLuwnusXg1xcVC3Lr/9BhcueGm4t2lj6kXr1zsOdepkPnTs2GFhu4RwIuEuvMeaNUXq7eCl4V6hgjl7qVi4g5RmhPeQcBfeodgiMomJULkytGhhcbtK07mzSXLb5PYmTaBqVQl34T0k3IV3cKq3gwn3jh3NdnZeqXNnM7k9KQmQQVXhfSTchXew74DdsSN5eWbijFeWZOyuv95crlvnONSpk9l2TwZVhTeQcBfeYe1axyIyu3aZTrFXh3uTJlCtWomDqjt3WtguIWwk3IX17PMebb1h+2CqfZDSKwUFmdJMsZ47SGlGeAcJd2G9YvMeExMhPByaNbO4XVfSuTNs22Z2jQKaNpVBVeE9JNyF9YrNe9y82VRogrz9f2fnzjgGCDDt7dBBwl14B2//8xGBwGneY0GBOdOzXTurG+UC+6Bqsbr7li2y/K+wnoS7sJ7TvMcDB+D8eR8J97p1zSqRTuGekADZ2TKoKqwn4S6slZdnNr6wlWS2bDGHfSLcwfTenQZV7dsBOu3lIYQlJNyFtZKSzICkU7gHBUHr1ha3y1WdO5sB4TNnADOoWqmSowwvhGUk3IW1ig2mbtliArJyZQvbVBb2urvt5wgOhrZtCz+BCGEVCXdhrcREM3+waVPAhKLPlGSg8Ewrp7p7+/am527by0MIS0i4C2slJpopJkFBnDsH+/ebnq/PiIoyZ6s61d3btYP0dDh82MJ2iYDnUrgrpW5WSu1WSiUrpSaWcPvDSqltSqnNSqkVSilv2xhNeKOcHNNVt/V+t241h32q5w6m/Rs3Oq62b28upe4urHTFcFdKBQPvA4OAVsDoEsJ7uta6jda6PfA68K9yb6nwPzt2mMVYfD3cO3Y0W+6dOgWYvTyUknAX1nKl5349kKy13qe1zgFmAMOc76C1Pud0tQog1UZxZSUMpkZFQb16FrbpahSb/xgeboYQZFBVWMmVcK8LOFcPU2zHilBK/UkptRfTc3+sfJon/FpiIkRGQuPGQOFgqlIWt6usOnQwl8VKM9JzF1ZyJdxL+lO7pGeutX5fax0P/BV4rsQHUupBpVSiUirx5MmTZWup8D+JiabXrhT5+WYNLp8ryQBER0PDhkXCvV072LcPzp61rlkisLkS7ilAnNP1ekDqZe4/A7itpBu01h9qrRO01gmxsbGut1L4n+xsk+a2kszevWYNd58MdzClGafTUu2DqvZxBCE8zZVwXw80VUo1UkpVBEYBs53voJRq6nR1MLCn/Joo/NL27WbLIlu92ueWHSiuY0fYswfOmeEnmTEjrHbFcNda5wGPAAuAJOBrrfUOpdTLSqmhtrs9opTaoZTaDDwJjHdbi4V/sPdyncI9OBha+eokWnvd3ZbmtWtDbKwMqgrrhLhyJ631PGBesWN/c/r3hHJul/B3mzZBRAQ0agSYEGzRAsLCLG7X1bLPmNm4EXr1QikZVBXWkjNUhTU2bSqyI8eWLT52ZmpxtWqZ7nqxGTP26pMQnibhLjwvP9+MNNpKGadPm1P1fbbebtex4yUzZi5ehN27LWyTCFgS7sLzfvvNTI2xhfu2beawX4R7UpL52SgcVJW6u7CChLvwPHsh2pZ+Pj9Txq5jRygocLxbNW8OoaFSdxfWkHAXnrdpE1Ss6Jgas2WLmVlSq5bF7bpWxc5UDQkx68xIuAsrSLgLz9u0yWy1VKEC4MPLDhRXv745W7WEZQhkbXfhaRLuwrO0NuFu6+Xm5ZkZJT5fkgHz7lTCoOqpU5B6uXO6hXADCXfhWSkpkJbmCPfffjMzSvwi3MGE+/btZq16ZFBVWEfCXXiW/cxUW7j7zWCqXceOJth37gQK5+5L3V14moS78KxNm0z5wpZ6W7aY0nuLFha3q7w4n6mKOQk3Pl7CXXiehLvwrE2boFkzs6MFJtxbtjSTZ/xCfLzZ8LtY3V3CXXiahLvwLKfBVCicKeM3goJMob1YuCcnQ0aGhe0SAUfCXXhOWprZa9QW7idPwtGjfhbuYEozW7aYZRYwP5/WhWfiCuEJEu7Cc+y1CVu4++yG2FfSsaNZguC334DCn0827hCeJOEuPKdYuPvdTBm7YoOqDRqYgVWZDik8ScJdeM6mTVCvHsTEACbs7Jta+BX7ojK2aZ/2yUES7sKTJNyF5/j7YKpdhQpmURmnPVXbtTNlGVmGQHiKhLvwjMxM2LXLEe7283z8MtzB/JxOi8q0awfnz8OBA9Y2SwQOCXfhGdu2meVwbefj79pldijy6d2XLqdDh8JdSCh8E5PSjPAUCXfhGf6+7EBx9kVlbD/3ddeZ2ruEu/AUCXfhGZs2QVSUmTqCCbnQUDP26JfatjVpbgv3KlWgaVOZDik8R8JdeIZ9Q2zbou1btpjebEiIxe1ylypVzDuX07oD7dpJz114joS7cL+8PFNzt5VktPbjmTLOOnQoMmOmbVvYu9cMrArhbhLuwv127YLsbEe4Hz9ulh7w+3Bv394st5CWBhT+vLIMgfAECXfhfoE2mGpnn9NvK83IMgTCkyTchftt3gxhYY7R04AJd/uMGVu4x8VBZKTU3YVnSLgL99u0yRScbaOnW7aYoIuKsrhd7hYbC3XryjIEwhIS7sK9tDY9V3svlgAZTLUrNqhqX4agoMDCNomAIOEu3OvQIUhPd9Sfs7PN+KrfnplaXIcO5gfOygJMuGdkwP79FrdL+D0Jd+FexZb53bnT7GERUD33ggLHFBlZhkB4ioS7cK9Nm8zWc23aAAE0mGpXwjIEQUES7sL9JNyFe23ebDbErlwZMKFWqRI0aWJxuzylYUMzRcYW7pUqmZdDwl24m4S7cK8SBlPbtIHgYAvb5ElKmZ+/2DIEMtdduJuEu3Cf9HQ4eDDwlh0orkMHk+ZOG2bv3w/nzlncLuHXJNyF+9h7q7ae+5EjJu8DLtzbtzezZXbvBgpnCknvXbiThLtwH3u429I84AZT7ezLENjq7rIMgfAECXfhPps3mx2wa9YECsPdNnEmcLRoYRavt73Z1a0L0dEyqCrcS8JduM/mzZdsiN2wIVSrZl2TLFFsw2xZhkB4goS7cI+LF80ZS4G67EBx7dubcHfaMHvbNscYqxDlTsJduMeOHWaTDlu4Z2bCnj0BHO4lbJidmQn79lncLuG3JNyFexSbKbN9uzkLP6DDHS5Z211KM8JdJNyFe2zaBOHhEB8PBPBMGbs2bYpsmN2qlTmRS8JduIvSthqgpyUkJOjExERLnlu4V25uLilffUV2nTpmtgymInHhglnH3bZHduBJTTVr2teoUdJVIQAIDg4mMjKSmJgYgoIu7X8rpTZorROu9Dj+uve8sFDK4cNUbdaMhjVroho0AMyqt9Wrm1mBASs01LzDtWxZ0lUh0FqTm5vL8ePHSUlJoX79+lf9WFKWEeUuOyOD6kFBKNtiYVqbEzQrVbK4YVarXBlycsxAM+b1cLoqBEopKlasSN26dcnIyLimx5JwF+UvJwcFjpUgc3LMlD/b1cBlf3fLzAQKXw/bPh5COJRUjinzY5RDO4QoKifHXNrCzB5e0nO3pbkt3ItlvRDlSsJdlL+cHJNctt6HPbx8KdwPHDiAUoo8W83k+PHj9OrVi6pVq/LUU09dcv977rmH55577vIPWqGC+bK9IBUqmAHV8uq5v/jii4wdO/aqvrdPnz58/PHH5dOQYl555RXuv//+Um+fMmUKPXr0cMtzBzIZUBXlLyenSA0mM9MMHvryGu4ffvghMTExnDt3DnUt030qV3akuVLmDc/fyzLPPvus498HDhygUaNG5ObmEhIi8eNO0nMX5evkSVNgd+qmZ2X5fr394MGDtGrV6tqCHQrDvaCgyFWLZiS7XZ6MFltGwl2UL/uZqbY0z883y8x4S0nm8OHD3HHHHcTGxlK9enUeeeQRAPLz83n66aeJiYmhcePGzJ071/E999xzD5999hmvv/464eHhLF68uMTHPnXqFAMHDqRq1ar07t2bgwcPOm6bMGECcXFxRLRoQadx41hue4xKlWDbtnV06pRAREQENWvW5Mknn3R835o1a+jWrRuRkZG0a9eOZcuWOW7bv38/vXv3pmrVqgwcOJBTp06V+nOnp6dz6623EhsbS1RUFLfeeispKSkl3jc/P5+nnnqKmJgYGjVqxL///e8iJarU1FSGDh1KdHQ0TZo04aOPPnJ874svvsidd97J2LFjiYiIYMqUKUXKRb169QIgMjKS8PBwVq9e7fjep59+mqioKBo1asT8+fMdx/v06cNzzz1Ht27dCA8PZ8iQIaSlpTFmzBgiIiLo3LkzBw4cKPVnD1haa0u+OnXqpIUfev11vXP+fK1zc7XWWp8/r/X69Vqnp1vcLq11Xl6ebtu2rX788cf1hQsXdFZWll6+fLnWWusPPvhAN2/eXB86dEinpaXpPn36aEDn2n6O8ePH60mTJpX62OPHj9fh4eH6l19+0dnZ2fqxxx7T3bt3d9z+xRdf6FOnTuncCxf0mxMm6JqxsTorK0tnZGjdpk0X/Z//fK611vr8+fN69erVWmutU1JSdHR0tJ47d67Oz8/XCxcu1NHR0frEiRNaa627dOmin3jiCZ2dna1/+eUXHR4erseMGVNi+06dOqW//fZbnZGRoc+dO6fvvPNOPWzYMMftvXv31h999JHjtWjZsqU+fPiwPn36tO7fv3+R16JXr176D3/4g87KytKbNm3SMTExevHixVprrV944QUdEhKiZ86cqfPz83VmZqZ+4YUXHO3av39/kcfSWuvJkyfrkJAQ/eGHH+q8vDz9n//8R9euXVsXFBQ42hYfH6+Tk5P1mTNndMuWLXXTpk31okWLdG5urh43bpy+5557XPo/4Et27txZ4nEgUbuQsVL0EuXLvmeqrZ4a9OTjNN+wmSpVcN/nxPbt4e23r3i3devWkZqayhtvvOGo99oH8r7++msef/xx4uLiAHjmmWeK9JJdMXjwYEfP9B//+AfVqlXj8OHDxMXFFQ50as1T48fz98mT2b17N23atCMkpAK//ZbMqVOniImJoUuXLgBMnTqVW265hVtuuQWAgQMHkpCQwLx58+jbty/r169n8eLFhIaG0qtXL4YMGVJq26pXr87vfvc7x/VJkybRt2/fEu/79ddfM2HCBOrVqwfAxIkTWbJkCWA++axYsYIff/yRsLAw2rdvz/33388XX3xB//79AejatSu33XYbAJVc/MjWoEEDHnjgAQDGjx/PH//4R44fP06tWrUAuPfee4m3LWUxaNAgdu7cyYABAwAYPnw4zz//vEvPE0ikLCPK16ZNULGi42peHqBAecH/tMOHD9OgQYMSB/JSU1MdwQ4mbMrK+fvDw8OJjo4mNTUVgLfeeouWLVtSLTKSyN69OXv+PKdOnSIoCP7+90/Ys+c3WrRoQefOnfnxxx8BU+f/5ptviIyMdHytWLGCo0ePkpqaSlRUFFWqVHGpzZmZmTz00EM0aNCAiIgIevXqxZkzZ8gvYc3h4q+F879TU1OJjo6matWqRZ73yJEjJd7fVfYQB6hsK+lduHDBcaymbcMXMG8Yxa8731cY0nMX5Scz0+wT6hTuR/78Nkp5x7IDcXFxHDp0iLy8vEsCvnbt2hy2LccLcOjQoTI/vvP3X7hwgdOnT1OnTh2WL1/Oa6+9xpIlS7juuusIOnKEqOuuQ9sGVVu0aMorr3xJ69YFfP/999x5552kpaURFxfHuHHjitS07Q4ePEh6ejoZGRmOgD906FCpA75vvfUWu3fvZu3atdSqVYvNmzfToUMHdAkjubVr1y5Sj3f+uerUqcPp06c5f/68I+APHTpE3bp1Hfe53KDzNQ9IC5d5QX9K+A37ur62cLcvO+AtM2Wuv/56ateuzcSJE8nIyCA7O5uVK1cCMGLECN59911SUlJIT0/n1VdfLfPjz5s3jxUrVpCTk8Pzzz/PDTfcQFxcHOfPnyckJITY2Fjy8vJ4+f33OZeR4TjZa968qRw/fpKCgiAiIyMBs3jU2LFjmTNnDgsWLCA/P5/s7GyWLVtGSkoKDRo0ICEhgRdeeIGcnBxWrFjBnDlzSm3b+fPnqVSpEpGRkZw+fZqXXnqp1PuOGDGCd955hyNHjnDmzBlee+01x21xcXF069aNZ555huzsbLZu3conn3zCmDFjXHqNYmNjCQoKYp8sZO92Eu6i/NiWs7WH+8WLJuu9ZaZMcHAwc+bMITk5mfr161OvXj2++uorAB544AFuuukm2rVrR8eOHbnjjjvK/Ph33XUXL730EtHR0WzYsIFp06YBcNNNNzFo0CCaNWtGgwYNCAsPJ65mTfMCAcuX/8TIkdcRGRnOhAkTmDFjBmFhYcTFxTFr1ixeeeUVYmNjiYuL44033qDA1uOfPn06a9euJTo6mpdeeom777671LY9/vjjZGVlOWr6N998c6n3feCBB7jxxhtp27YtHTp04JZbbiEkJIRg24kKX375JQcOHKBOnTrcfvvtvPTSSwwcONCl16hy5cpMmjSJ7t27ExkZyZo1a1z6PlF2suSvKD9/+AN8+SVJq1fTsmVL0tNh716z6qFTaVgUFJg3who1IC6O3FyzrntcnGMvca8yf/58Hn744SJTO4X7JSUl0bKEJUNdXfJXeu6i/Nhnytj44rIDHhEUVOTU1PJehuBaZWVlMW/ePPLy8jhy5AgvvfQSt99+u9XNEmUk4S7KR34+bN16SbiHhTmWmBHOKlc2L5Dtk7P9qjfQWvPCCy8QFRVFhw4daNmyJS+//LLVzRJlJLNlRPnYs8ekk32vUExPNDzcwjZ5s8qV4dQpyM2FihWpVAlOnDBZb/WEksqVK7N+/XprGyGumfSpRPkotiF2bu4l64cJZ8XIsotxAAAgAElEQVSW/61c2QR7draFbRJ+RcJdlI/Nm03x2DYAZK8fS7iXothi7var3lJ3F77PpXBXSt2slNqtlEpWSk0s4fYnlVI7lVJblVJLlFJlP71P+LZNm+C66xzTIGUw9QqCg82AhO2FCgsz5RhvqbsL33fFcFdKBQPvA4OAVsBopVSrYnfbBCRordsC3wKvl3dDhRfTGjZuhI4dHYcyMwv3phClqFTJkeZBQSbgpecuyosrPffrgWSt9T6tdQ4wAxjmfAet9VKttb3PsQaoV77NFF7t8GEzONipk+OQN52Z6rWKbZjtTTNmhO9zJdzrAoedrqfYjpXmPmD+ZW4X/mbjRnNp67l727IDXquEPVVzc80XXNu2ed5GKUVycnK5P65L2xtepWnTpnHjjTeWevuyZcscK2d6I1fCvaSJWSWe1qqUGgskAG+UcvuDSqlEpVTiyZMnXW+l8G4bNpgacrt2QOH+2BLuV1DCjBmQ0oy3GDNmDAsXLnRcd9cblLu4Eu4pgPManvWA1OJ3UkoNACYBQ7XWF0t6IK31h1rrBK11Qmxs7NW0V3ijjRvNLBnb6Kk93ANpMPWqtpMrtmG2t8yYka3x/OM1cCXc1wNNlVKNlFIVgVHAbOc7KKU6AP/DBPuJ8m+m8Fpam567U709N9cMEIaGWtiuUjRs2JA333yTtm3bUq1aNUaOHEm2bXK5/WP2K6+8QkxMDA0bNnQs/lXaY7322mu0bduWKlWqkJeXx6uvvkp8fDxVq1alVatWzJw503H/KVOm0KNHj6Lbya1f70jzlJT9PPxwb+LjS942b/bs2Vx33XVERkbSp08fkpKSirTljTfecLTlvvvu4/jx4wwaNIiqVasyYMAA0tPTS/w5pkyZQvfu3XniiSeIjo7mxRdfZO/evfTr14/q1asTExPDmDFjOHPmjEuvI8Abb7xB7dq1qVOnDp9++mmR57vnnnv405/+xODBg6latSo33HADe/fuLfV1Hj58OLVq1aJatWr06tWLHTt2lHrf119/3fG8H3/8cZHe9tmzZ7n77ruJjY2lQYMG/P3vf3cswlbSa2D/fUHh9oDt2rUjPDzcseAcmOWUa9SoQe3atZk8eXKRn/OPf/wjgwYNIjw8nO7du3Ps2DEef/xxoqKiaNGiBZvsi+25gyvbNQG3AL8Be4FJtmMvY8IcYDFwHNhs+5p9pceUbfb8xJEjWoPW77zjOPTzzzt1UpKFbbqMBg0a6M6dO+sjR47otLQ03aJFC/3BBx9orbVeunSpDg4Odmxdt2zZMl25cmW9a9euUh+rXbt2+tChQzozM1NrrfXXX3+tjxw5ovPz8/WMGTN05cqVdWpqqta6lO3katbUBevWaZ2fr7t06aLvuecJvWnTpdvm7d69W1euXFkvXLhQ5+Tk6Ndee03Hx8frixcvOtpyww036GPHjumUlBQdGxurO3TooDdu3Kizs7N137599YsvvljizzF58mQdHBys3333XZ2bm6szMzP1nj179MKFC3V2drY+ceKE7tmzp54wYYJLr+P8+fN1jRo19LZt2/SFCxf06NGjNaD37NmjtTZbEkZFRem1a9fq3Nxcfdddd+mRI0eW+jv75JNP9Llz53R2draeMGGCbteuneM25+0P58+fr2vWrKm3b9+uMzIy9NixY4s877hx4/TQoUP1uXPn9P79+3XTpk31xx9/XOprMHny5CJbJTo/lvP/l+eff17n5OTouXPn6kqVKunTp0872la9enWdmJios7KydN++fXXDhg31Z599pvPy8vSkSZN0nz59Sv25r3WbPdlDVVyb2bPNf6MVK7TWWufna/3TTzv1wYPm5gkTtO7d271fTplzRQ0aNNBffPGF4/qf//xn/dBDD2mtC/9YL1y44Lh9+PDh+uWXXy71sT755JPLPl+7du30Dz/8oLU2ARIfH++4LSMjQwP66Pz5+mBSkg4ODta7d1/QiYnmdRw9erQj3F9++WU9fPhwx/fm5+frOnXq6KVLlzraMnXqVMftd9xxh3744Ycd1999990ie6Y6mzx5so6Li7vszzFz5kzdvn37Ij97aa/jvffeq//61786btu9e/cl4X7fffc5bp87d65u3rz5ZZ/fLj09XQP6zJkzjseyh/u9996rJ06c6Ljvnj17HM+bl5enK1asqHfs2OG4/b///a/u3bt3qa+BK+EeFhZWZD/Y2NhYxx6448eP1/fff7/jtnfffVe3aNHCcX3r1q26WrVqpf6s1xrucoaquDYbN5qzb2yDqXv3mkqNNw+mFt/SzXmLtpK2rrNvlVeS4lvKff7557Rv396xLd727duLlFdK3E4uK4vU/fuJioqievUqaG2WenfeNi81NbXI9aCgIOLi4opsb3ctW9EV/zlOnDjBqFGjqFu3LhEREYwdO/aSMlFpr6MrWxZe7nfgLD8/n4kTJxIfH09ERAQNGzYEuKQtJT2v879PnTpFTk5OkbaUx/aA1atXL7KrV/GfxcrtAWXhMHFtNmwwe+jZVgjbvBmqVi0cHHRh32qvUtLWda1bty71/s7bxh08eJAHHniAJUuW0LVrV4KDg2nfvn2JW9kVERRE7YgI0tPT0ToDqEJmZtFt8+rUqcO2bdsc36K15vDhw0W2t7sWxbe/e+aZZ1BKsXXrVqpXr84PP/zAI4884tJjlceWhXbTp09n1qxZLF68mIYNG3L27FmioqLKvD1gTEwMFSpU4ODBg7Rq1crRLle3B/RF0nMX12bjxiKDqfb1w3x5pox967rly5fz448/Mnz4cJe+LyMjA6UU9plgkydPZvv27Vf+xrAwGkRHk5CQwD//+QJ5eTn88kvRbfNGjBjB3LlzWbJkCbm5ubz11luEhobSrVu3q/oZr+T8+fOEh4cTGRnJkSNHeOONEmc3l2jEiBFMmTKFnTt3kpmZedkt/VxpR2hoKNWrVyczM5Nnn332ss87efJkkpKSyMzMLLJMcXBwMCNGjGDSpEmcP3+egwcP8q9//atM5xHUrFnTp7YHlHAXV+/4cThypMiyA/b1w3x1DfdatWoRFRVFnTp1GDNmDP/9739pYdvd+5VXXmHQoEGlfm+rVq146qmn6Nq1KzVr1mTbtm107979yk9q27hj+rRprFu3ln79onnrraLb5jVv3pypU6fy6KOPEhMTw5w5c5gzZw4VnTYjL4vw8HCWL19e6u0vvPACGzdupFq1agwePLhM2w4OGjSIxx9/nH79+tGkSRP69et3VW0EuPvuu2nQoAF169alVatWdOnS5bLP+9hjj9G3b1+aNGlC165dAQi1Tdt67733qFKlCo0bN6ZHjx7cdddd/P73v3e5LS+++CLjx48nMjKSr7/++qp/Jk+RbfbE1Zs3DwYPhl9+AdtUsTp14KuvkujZ89LtwbzdsmXLGDt2bJGP9h5x6hQcOGAWXqtUif374dw5xzCGuEpJSUm0bt2aixcvFqmL+wrZZk9Yx77sgG0N9+PH4ehRx8KQwlUlnKnqvAyBcN3MmTPJyckhPT2dv/71rwwZMsQng708SLiLq7dhAzRrBhERgNnkGSTcy8y+3q/tZKZiWS/K4H//+x+xsbHEx8cTHBzMBx98YHWTLBOYb2mifGzcCE4DevbBVF8N9z59+ni+JAOFG2YXW4YgMxOqVfN8c3zZTz/9ZHUTvIb03MXVOXUKDh26ZKZM/fq+O5hqKacNs0NCzNIN0nMX10L+DMXVKbbML5jNmOz7Y1s1UO+zKlc267rbVl2Ttd0DW3n8/Ui4i6uzYYO5tIV7Rgb89puZ4REcHEyujAaWjf2s2IwMwIT7xYuOfTxEgMnKyqLCNW5jJuEurs7GjRAfD5GRgCnJFBRAQgJERkZy/Phxx4p7wgWVKhXZRNWe9dJ7DyxaazIzMzly5Ag1atS4pseSAVVxdTZsMEnudBVMCT4mJoaUlBR2795tUeN81NmzcP48nD9Pfr4Z1sjLk0HVQFOhQgVq1qxJhG0W2tWScBdll54O+/fDQw85Dm3YALVqmZOYIIj69etb1jyf9d57MG2aeX2Dgrj1VujSBWbMsLphwhdJWUaUXQmDqcX26xBXo3Nnc2rqnj2AeXntL7UQZSXhLspu/XpzaUvzjAxISpJwv2b2Mpft9e3UyeT8uXMWtkn4LAl3UXbr10OTJhAdDRQOpkq4X6OWLc00GduaS/YPRvaTw4QoCwl3UXbr1sH11zuuOg+mimsQEmIS3dZzt4e7/fUVoiwk3EXZpKZCSsol4V6zpn0wVVyThARzNlheHjVrQt26UncXV0fCXZSNvd7eubPjkH1WpJ9tZGONzp3NAmI7dwIyqCqunoS7KJt16yA42LHOgAymljP7m6ZTaWbXLseJq0K4TMJdlM369dC2rWPpwi1bZDC1XMXHm7OWbIOqnTqZ13fTJovbJXyOhLtwXUGBCXcZTHWfoCBT47L13It15IVwmYS7cF1yMpw5I4Op7paQAFu3wsWL1KoF9epJuIuyk3AXrlu3zlwWG0zt1EkGU8tV585mj72tWwHzXirhLspKwl24bt06s1xhq1aAWbFw504pyZQ7+5vn2rWOq8nJcPq0hW0SPkfCXbhu/XqT5MHBgOm1FxQUqdKI8hAXB7VrFwl3cIyxCuESCXfhmpwcM2XDKcntVRoJ93KmlFkOcs0aoPCTkZRmRFlIuAvXbNtmtgYqFu4NG8I17ikgStK1q6nFnDpFZCQ0b174ZiqEKyTchWtKGExdu1Z67W7TpYu5tPXeO3c2vwLZmla4SsJduGbdOoiNhQYNADh+HA4elHB3G/vYhlO4HzsGR45Y3C7hMyTchWvs3XTbnEd7/feGGyxskz+rXNnsNm4Ld/ubqNTdhask3MWVnT5tFpDp1s1xaO3aIkvMCHfo0sV8YsrPp317syKwhLtwlYS7uDJb79E53Netg9atzbR34SZdu5oNs5OSCAuDNm1kUFW4TsJdXNmqVaabbhtMLSgwISMlGTezD6quXg2Ylz8x0bz+QlyJhLu4stWroX17Rze9hCVmhDvEx0P16kXq7mfPOvbPFuKyJNzF5eXlmQJ7166OQ/bSgPTc3azYyUz219teJRPiciTcxeVt22Z2iig2mFqlitnPWbhZly5mAZ8zZ2jVCiIiHFUaIS5Lwl1c3qpV5rLYYGpCgmOJGeFO9rr7unUEBZneu4S7cIWEu7i8VavMYu316wNmBYLNm6Uk4zH2cwtstZhu3cyHqXPnLG6X8HoS7uLyVq0yiWI7eWnjRrOGmIS7h0REwHXXObrrXbuaJQhkSqS4Egl3UbqjR+HAgSKDqSVUaYS7de1qwj0/3/GmKqUZcSUS7qJ09gRxSvKVK6FxY6hVy6I2BaKePc0cyO3biYws0pEXolQS7qJ0q1ZBaKhjjQGtTbh3725xuwJNr17mcvlyoLAjLyczicuRcBelW7XKTIsJDQVg7144cULC3eMaNDC7MzmF+5kzsHu3xe0SXk3CXZQsO9vso+dUkrHX2yXcLdCzpwl3rR2/EinNiMuRcBcl27DBTItxGkxduRKqVXPsjy08qWdPM8C9bx/NmkFUVOGbrRAlkXAXJfvlF3PZs6fj0MqVJuuD5H+N59l/D8uXExRUWHcXojTyZypKtmyZWdM3JgaA9HTYsUNKMpZp2RKio4vU3W2rEghRIgl3cancXPOZv3dvx6ESZkUKTwoKgh49ioQ7SO9dlE7CXVxqwwazWJhTuNuXdJczUy3Us6dZ7/fYMbp0MTsz2bJeiEtIuItL2evtTuG+cmWRJd2FFZzq7lWqmFmq9l+VEMVJuItLLVtmarw1agCmSrN2rdTbLdexo9k429Zd79XL7KmamWlxu4RXknAXReXlwYoVRXrtmzZBVpbU2y1XoYIpttvCvXdv88Yrm3eIkki4i6I2bYILF4qE+7Jl5tLpkLBKz56wZQucPUv37mac9ddfrW6U8EYS7qKoEurt9iqNLBbmBXr1Mov8rFhBtWpmHETq7qIkEu6iqGXLoFkzqF0bMFWa5cuhTx9LWyXsunQxa/0sXQqYrF+zxmyiIoQzCXdRKD/fJLlTr33DBlOlkXD3EpUqmcGPJUsA86vKzobERIvbJbyOhLsotGWL2b9N6u3erV8/s9dhWho9ephDUpoRxUm4i0KXqbfXrGlNk0QJ+vUzl8uWERNjNu+QQVVRnIS7KLR0KcTHQ716gJlmt2KFlGS8TufOEB4OP/8MmPfilSvN+IgQdhLuwsjNNd30AQMchzZuNPX2vn2ta5YoQYUKZiTVVnfv1cv8njZutLhdwqtIuAtj7Vo4fx4GDnQcsk3IkHq7N+rXz2zFdOSI4/dj/30JAS6Gu1LqZqXUbqVUslJqYgm391JKbVRK5Sml7iz/Zgq3W7TInBFjr+diOvKtWjlWIRDexP57WrqUWrWgTRvzKxTC7orhrpQKBt4HBgGtgNFKqeJ78RwC7gGml3cDhYcsXmxWooqKAqTe7vXatYPq1R2JPmCA+X1lZVncLuE1XOm5Xw8ka633aa1zgBnAMOc7aK0PaK23ArIfuy86e9aUZZzq7YmJZtVfCXcvFRRkfl8LF4LWDBxoTmSSJYCFnSvhXhc47HQ9xXZM+Itly8wJTE719gULQKkiVRrhbW66CY4dg61b6dULKlaU0owo5Eq4qxKO6at5MqXUg0qpRKVU4smTJ6/mIYQ7LFpklpJ12gx7wQIz4656dQvbJS7vxhvN5YIFVKliTlyVcBd2roR7ChDndL0ekHo1T6a1/lBrnaC1ToiNjb2ahxDusGiRmRITGgqY/VLXrTMdQ+HF6tY1I6kLFgDmg9eWLXD8uMXtEl7BlXBfDzRVSjVSSlUERgGz3dss4TH79sFvvxX2AjFjqwUFEu4+4aabzEhqRoajqmab/i4C3BXDXWudBzwCLACSgK+11juUUi8rpYYCKKU6K6VSgOHA/5RSO9zZaFGO5s0zl4MHOw4tWADVqsl+qT7hppsgJweWLaNjRzPZafFiqxslvEGIK3fSWs8D5hU79jenf6/HlGuEr5k3D5o2NV+YpcIXLID+/c0GzMLL9ehhVopcsIDgwYPp399U2bQ2A+IicMkZqoEsM9Oc1njLLY5DSUmQkgI332xhu4TrwsLM+hDz5zumRKakmJNXRWCTcA9kS5eaxcCdwt02Nif1dl8yeDAkJ8Pu3Y66u/33KAKXhHsgmzfPTIF0WjxmwQJo0QLq17ewXaJsbr3VXM6ZQ6NGZonmOXOsbZKwnoR7oNIa5s41ZznapkBmZZkl3aXX7mPq1zfLEdgSfcgQ83s8d87idglLSbgHqp074eDBIiWZX381VRoJdx80ZIhZ1D0tjSFDzNruUpoJbBLugWrWLHNp/0hvO1Sliqwn45OGDDEnJ8yfT9euEB0tpZlAJ+EeqGbOhOuvN2c5YnJh1izTa69UyeK2ibJLSIBatWDOHIKDzQeyefPMkkEiMEm4B6LDh82yj7ff7ji0YQOkpsJtt1nYLnH1goLMrJmffoKcHIYMgbQ0WL3a6oYJq0i4ByJ7ScYp3H/4AYKDi5yoKnzNsGFmFHXZMm66yZyEJqWZwCXhHohmzjTzHZs3dxyaNcvsxRkdbWG7xLUZONBsnP3tt1SrZma4SrgHLgn3QJOWZubJOfXak5Nhxw4pyfi8sDAzsDpzJuTlMWSIOeN4716rGyasIOEeaH780YyyOYW7vUozbFgp3yN8x513wqlT8OuvjolQs2UN14Ak4R5ovvkG4uKgUyfHoR9+gPbtoUEDC9slysfNN5uzjr/9lvh4c27TN99Y3ShhBQn3QHL6tDmzZeRIM7sCOHHCnPsiJRk/UbmyGRX//nvIz2fkSDNj5tAhqxsmPE3CPZB8/705dXHUKMehWbPMSgRSkvEjv/ud2Y5p5UpGjjSHvv7a2iYJz5NwDyQzZkCTJtCxo+PQ9Olm0ky7dha2S5SvwYPNmWgzZtC4sTm/6auvrG6U8DQJ90Bx7JhZ4nfUKMcuDikpZuLMmDGysYNfCQ83dbavvoKcHEaMMOesyayZwCLhHii+/dasMTB6tOPQjBmmJON0SPiLsWPNGMtPPzFihDkkpZnAIuEeKL78Elq3hlatHIemTzfLyzRpYmG7hHsMHAixsTB1Kg0aQJcuUpoJNBLugeC332DVKhg3znEoKQk2bYK77rKwXcJ9KlQwJbjZs+HsWUaOhC1bZPu9QCLhHgimTDELxziF+/TpZjakfTaF8ENjx8LFi/DddwwfbsZVZsywulHCUyTc/V1+Pnz+uTm5pXZtwNTZp0+H/v3NKrHCT3XuDE2bwmefUbcu9Otn3ucLCqxumPAECXd/t2gRHDkC997rOLRuHezbJyUZv6cU/P73Zout3bu5/344cAB+/tnqhglPkHD3d5MnQ/XqZkEpm08+MdOgnZaXEf7qnnvM2r8ff8xtt0FUFHz8sdWNEp4g4e7P0tLMwjFjxkDFioBZ7nv6dDP9sVo1i9sn3K9WLRg6FKZMIUxdZOxYs2hkWprVDRPuJuHuzz79FHJy4P77HYemT4eMDHjoIQvbJTzrgQfMSpGzZnHffea/xLRpVjdKuJvSWlvyxAkJCToxMdGS5w4I+flmMK1+fVi2DDADqR06mFLsxo1yVmrAyM+Hxo2hWTNYtIjOnc0kmi1b5P+AL1JKbdBaJ1zpftJz91c//QT798Of/uQ4tG6d+YN++GH5ow4owcGm9754MezaxX33wbZtZkkC4b8k3P3V++9DnTpF1vL973/NsiMySyYAPfgghIbCu+8yerQZUP/f/6xulHAnCXd/lJwM8+ebP+gKFQA4c8acfn7XXVC1qsXtE55Xo4YZWP/sM6rln+buu2HqVLOev/BPEu7+6O23Tag/+KDj0JQpkJUlA6kBbcIEyMyEjz7i8cdN3f0//7G6UcJdZEDV35w8afbLGz3aTGgHcnMhPt4cXr7c4vYJa/Xvb9Ya2rePIXdUYO1aOHjQlGmEb5AB1UD13nuQnQ1//rPj0JdfwuHD8MwzFrZLeIcnnjAL+X/1FU89ZfoCU6da3SjhDtJz9ycXLpipj336mC31MOuItGljJkzI1DdBQYHZDT03F71tO52uDyY7G7Zvd2yrK7yc9NwD0UcfQXo6/PWvjkM//gg7d8LEiRLsApPgzz0Hu3ahvvuWp54yyz//9JPVDRPlTXru/iIjwxTWW7VyrAylNXTvDkePwp49ZokRIcjPd3ycy03cQuMmQdSvDytWSAfAF0jPPdD8+99mx/u//91xaPlyWL0ann5agl04CQ6GSZNg+3YqzP2B554ze7nMn291w0R5kp67Pzhzxpxe3q2bqcNgeu39+pmSzP79ULmyxW0U3iUvz2y7GBRE7oattGgdQkQEbNggtXdvJz33QPKvf5lau1OvfdEis6TMc89JsIsShITAP/8JSUlU+OJTXn4ZNm82+6gL/yA9d1+XmmoWhLrlFsf29gUFZuPrU6fMnpmhoRa3UXgnraFnT0hOJn93Mu17hJOTAzt2SBnPm0nPPVD85S/mI/arrzoOTZ9uPl6//LIEu7gMpeDNN+H4cYL/703+/ndzftNnn1ndMFEepOfuy1asMD2vSZMcJZmMDGje3KwZtmaN1E+FC0aMgB9/RG/fQbcxjdi/H3btgshIqxsmSiI9d3+Xnw+PPgr16hU59fTVV82WqW+/LcEuXPTWWxAcjHr0Ed7/t+bkSdNfEL5N/vx91bvvmhGwN9+EKlUA09t67TUYO9ZMnBHCJXFxpoY3bx4d93/Ho4/CBx+Y9f+F75KyjC/aswfatoUBA2D2bFAKraFvX9i61YR8jRpWN1L4lLw86NwZjh/n3JqdtOwaSY0asH69DK56GynL+Kv8fLj3XggLM7st2E4p/O9/4Zdf4PXXJdjFVQgJMctXnDhBxLOP8M475oPhe+9Z3TBxtSTcfc3//R+sXAnvvGNGTYG9e81ZqDfeCPfdZ3H7hO9KSIC//Q2mTeN3OV8yeDA8+6zZkk/4HinL+JLVq6FXLxgyBL77DpQiN9ccSkoyK/vVq2d1I4VPy8sz/6F27uT44m20HxJHVJQpz9iGdoTFpCzjb06dMlPW4uLg008d5ZhJk8yUx48/lmAX5SAkxCzwXlBAzT/cwdRPLrJrFzz2mNUNE2Ul4e4L8vLMFJgTJ+CbbxwTkL//Ht54A/74R7jzTovbKPxH48Ym4BMT6f/Vgzz7jObTT2HaNKsbJspCwt3baQ2PPAILFsD770OnToCZFXP33dCli5mmLES5GjoUXnoJPv+cFyPfpkcPeOABUxkUvkHC3du9/rqZFTNxItx/P2B2SRs8GKpVM733sDCL2yj803PPwe9+R8hfnuS70d9Sp44Z7tm92+qGCVdIuHuzDz80oT5qFPzjHwCkpcHNN8PZszB3LtSubXEbhf8KCjLlmT59qDFhNAueWUZQkPn/d+yY1Y0TVyLh7q0+/BAeesis9jhlCgQFceyY2R41ORl++MFshSmEW4WFwaxZ0LYt8Y8MYt6zKzh5Em66yewNI7yXhLu30drMZbcH+/ffQ2goKSnQuzfs22d67P36Wd1QETAiIsyYz3XXkfDnvsx8bCnJydCjBxw4YHXjRGkk3L1JXp5ZDOzJJ+GOOxzBvmuXmXp89Kj5G+vf3+qGioATE2P25u3WjYGv9mfxvdNIS9N0727OrxDeR8LdW5w4YUZJ33/fnG76zTcQGsp335klPy5cgCVLTG9JCEtERMBPP8Hw4XR9fyzLu/4FdAHdu8sOTt5Iwt0b/PwztGsHv/5q1vd44w1y84P4y1/M/PXrroONG03IC2GpSpVgxgz45z+5bv5brK4ykBZxGQwfbj50XrxodQOFnYS7ldLTTW29f39zYtK6dXD//axbZ5b5sJ+g9Msvcvap8CJKmVlcCxdSP2s3y3fF8mSXVfz73+a8C1kq2DtIuFshL8+sF9CypWiTKvEAAAfWSURBVLl86ilITCS9XhsefdT8gaSlwcyZpkojW+UJrzRgAGzbRsVRd/DWmu7Mrnk/xw9m0aWL5sEHzYoZwjoS7p6Umwtffglt2pjT/Ro1gvXrSZ/0Ji+8XoWGDU2Y/+lPsHMn3Hab1Q0W4gqiosxc+CVLGBKzhl3pNXmyxjQ+/aSApk01L7xgOirC8yTcPeH4cVNjadIE7rrLnBwycyZJn6ziyakdadTIbIQzcGDhGtoREVY3Wogy6NcPNm8m4qN/8Wal59lc0Ja++Yt5+WVo0EDzxBOmwyI8R5b8dZeTJ82E9O+/h/nzHUupHrn3OWZn9Gfal0GsXAkVKphZj88+azZXEsLn5ebC9Onw3nvs2JDFa0HPMl2PIl8Hc31CPuPvDWbYMKhb1+qG+iZXl/x1KdyVUjcD7wDBwMda61eL3R4KfA50AtKAkVrrA5d7TL8L93PnYO1as5HGokVmhSWtOVe7Oau6P83y6NtYuDEG+4/cvLnZWGP8eNk5Sfgprc3o6rRpHJ+xlGknBzJZ/Z7tujUA7Zte4Jbbw+jVL4QbbnAsdiquoNzCXSkVDPwGDARSgPXAaK31Tqf7/BFoq7V+WCk1Crhdaz3yco/rs+F+7hwcOmS2P9q+HXbsgG3byN6ezCHi2KHasL32ALZH9WR7Vjy7DoRRUKAIDjZTGYcOhWHDzFiqbUl2Ifxffj6sXIn+fiY75h9i7m9NmMtgVtGNfEJQFNCq1mnaNM2mZduKtOpajZbtQ2nSRCYUFFee4d4VeFFrfZPt+jMAWut/Ot1nge0+q5VSIcAxIFZf5sG9Ldzz8iD713Vk7z9K9olzZJ88z8VT58lKy+RcWi7px7I5fSKf9KxQ0okinSiOUpvDFeM5TByncqoVebzGjaF1azN9vWdP6NoVwsMt+uGE8DZpafDrr5xfsYV1yy+yamckazJas5NWHKBRkbvGxpoSjv0rJsb08ot/hYebpXCcv0JDzRCXP3E13F3Z17wucNjpegpwQ2n30VrnKaXOAtWBcp8M9cUX8O67UFBQ+pfWl7/dfp/8fMjJgexs82+43qU2VAjOJyqigJp1gohrEMz1cWYeev36pkfesqUEuRCXVb063H47VW+/nf5AfzBzJ/ftI2NnIrvXnyNpRz57K7XhSP2uHDkCR46Y7f7S0ux/r66pWNGEfHCw+QoKuvK/g4OLfrK2/7ukY1dz+1/+Ysba3MmVcC+peFC8R+7KfVBKPQg8CFC/fn0XnvpSlSqZGnVQUMlfSpV+W0m3h4Y6vdOnHyUsVBMWXZmw6pUJq1qR0FAzcyUqCqKjzWXlysGYapUQotzExEBMDFWuv56O90DHUu6mNWRkwJkz5is93VxmZpqOWmlf+fmmY5eff+m/S7ru/HzOl6X9uyy3e6LU5Eq4pwBxTtfrAaml3CfFVpapBpwu/kBa6w+BD8GUZa6mwXfe6c4t5WRxdCG8nVLmk3F4uJy5fTmuVKPWA02VUo2UUhWBUcDsYveZDYy3/ftO4OfL1duFEEK41xV77rYa+iPAAsxUyE+11juUUi8DiVrr2cAnwBdKqWRMj32UOxsthBDi8lwpy6C1ngfMK3bsb07/zgaGl2/ThBBCXC0/myQkhBACJNyFEMIvSbgLIYQfknAXQgg/JOEuhBB+SMJdCCH8kIS7EEL4IQl3IYTwQxLuQgjhhyTchRDCD1m2h6pS6iRw0I1PEYMb1pP3QfI6FJLXopC8FoV87bVooLWOvdKdLAt3d1NKJbqyW4m/k9ehkLwWheS1KOSvr4WUZYQQwg9JuAshhB/y53D/0OoGeAl5HQrJa1FIXotCfvla+G3NXQghApk/99yFECJg+X24K6WeVkpppVSM1W2xilLqDaXULqXUVqXUTKVUpNVt8jSl1M1Kqd1KqWSl1ESr22MVpVScUmqpUipJKbVDKTXB6jZZSSkVrJTapJT60eq2lDe/DnelVBwwEDhkdVsstghorbVuC/wGPGNxezxKKRUMvA8MAloBo5VSraxtlWXygKe01i2BLsCfAvi1AJgAJFndCHfw63AH/g/4CxDQAwta64Va6zzb1TVAPSvbY4HrgWSt9T6tdQ4wAxhmcZssobU+qrXeaPv3eUyw1bW2VdZQStUDBgMfW90Wd/DbcFdKDQWOaK23WN0WL/N7YL7VjfCwusBhp+spBGigOVNKNQQ6AGutbYll3sZ0/gqsbog7hFjdgGuhlFoM1CrhpknAs8CNnm2RdS73WmitZ9nuMwnzsXyaJ9vmBVQJxwL605xSKhz4Dnhca33O6vZ4mlLqVuCE1nqDUqqP1e1xB58Od631gJKOK6XaAI2ALUopMGWIjUqp67XWxzzYRI8p7bWwU0qNB24F+uvAm/+aAsQ5Xa8HpFrUFssppSpggn2a1vp7q9tjke7AUKXULUAYEKGUmqq1Hmtxu8pNQMxzV0odABK01r60OFC5UUrdDPwL6K21Pml1ezxNKRWCGUjuDxwB1gN3aa13WNowCyjT2/kMOK21ftzq9ngDW8/9aa31rVa3pTz5bc1dFPFvoCqwSCm1WSn1X6sb5Em2weRHgAWYAcSvAzHYbboD44B+tv8Lm229V+FnAqLnLoQQgUZ67kII4Yck3IUQwg9JuAshhB+ScBdCCD8k4S6EEH5Iwl0IIfyQhLsQQvghCXchhPBD/x9TfqA9Z8GLLwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "import seaborn\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "\n",
    "n_samples=200\n",
    "np.random.seed(42)\n",
    "start=time.time()\n",
    "uniform_samples=np.random.rand(n_samples)\n",
    "gaussian_samples=np.sqrt(2)*scipy.special.erfinv(2*uniform_samples-1)\n",
    "end=time.time()\n",
    "print(\"The time of sampling %d samples using inverse cdf algorithm is %.6f\"%(n_samples,end-start))\n",
    "\n",
    "start=time.time()\n",
    "gaussian_samples_randn=np.random.randn(n_samples)\n",
    "end=time.time()\n",
    "print(\"The time of sampling %d samples using np.random.randn is %.6f\"%(n_samples,end-start))\n",
    "plt.figure(figsize=(6,6))\n",
    "seaborn.distplot(gaussian_samples,\n",
    "                 hist=False,\n",
    "                 kde=False,\n",
    "                 fit=scipy.stats.norm,\n",
    "                 fit_kws={\"color\":'r'},\n",
    "                 label=\"cdf based algorithm\"\n",
    "                )\n",
    "seaborn.distplot(gaussian_samples_randn,\n",
    "                hist=False,\n",
    "                kde=False,\n",
    "                fit=scipy.stats.norm,\n",
    "                fit_kws={\"color\":'b'},\n",
    "                label=\"np.random.randn algorithm\")\n",
    "plt.legend(fontsize=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 Sampling from a Gaussian (Box-Muller method)\n",
    "The inverse cdf algorithm is pretty slow for Gaussian distribution. We can use the Box-Muller method for faster sample. Suppose $U_1$ and $U_2$ are independent samples chosen from the uniform distribution on the unit interval $(0,1)$. Let $Z_0=\\sqrt{-2lnU_1}cos(2\\pi U_2)$ and $Z_1=\\sqrt{-2lnU_1}sin(2\\pi U_2)$.Then $Z_0$ and$Z_1$ are independent random variables with a standard normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The time of sampling 200 samples using Box-Muller algorithm is 0.000191\n",
      "The time of sampling 200 samples using np.random.randn is 0.000125\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAFpCAYAAABnHGgVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4lFX2wPHvTSONJNRQpbdICCWBIEgoglIFBQFdRUGsqLjoLr/VFda2q6KyuHZdWBVEehGkiBSpSRBCr9JFCB2SkHp/f9xkzJCQApO8M5PzeR6eYd55886ZTHJy59ymtNYIIYRwLx5WByCEEMLxJLkLIYQbkuQuhBBuSJK7EEK4IUnuQgjhhiS5CyGEG5LkLoQQbkiSuxBCuCFJ7kII4YYkuQshhBvysuqJK1eurOvWrWvV0wshhEvavHnzGa11lcLOsyy5161bl/j4eKueXgghXJJS6khRzpOyjBBCuCFJ7kII4YYkuQshhBuS5C6EEG5IkrsQQrghSe5CCOGGJLkLIYQbkuQuhBBuSJK7EEK4IUnuQgjhhiS5CyGEG7JsbRkhrPDrr3D0KNSoAY0agVJWRyREyZCWuygTtmyB6Gho0AC6dIEmTaBjR1i/3urIhCgZktyF25s2Ddq3h+PHYcIEWL4cJk6Ew4chJgZmzrQ6QiEcT8oywq399BM89JBppc+cCVWyV8G+4w54+GHo0weGDAEfH7j7bktDFcKhpOUu3NbhwzBwIDRtCgsW/JHYcwQHw5Il0KoVPPIIHDtmSZhClAhJ7sItaQ1PPw1paSaxBwXlf15AAHz7LaSnw4MPmq8Twh1Ichduae5cWLwYXnsN6tcv+NxGjeDdd2H1aqm/C/dRpOSulLpLKbVXKXVAKTW2gPMGKqW0UirScSEKUTzp6TBmDLRoAc88U7SvGTECwsNh7FhITS3Z+IQoDYUmd6WUJ/Ah0BMIA4YqpcLyOa888CywydFBClEcX39t6u3//Cd4FXHIgKenGUlz6BB8+mmJhidEqShKy70tcEBr/avWOg2YDuQ3ruA14G3gqgPjE6JY0tPhjTcgMhJ69ize1/boAbffbko0GRklE58QpaUoyb0mkHscwfHsYzZKqVZAba319w6MTYhimznTzEL9+99vbPbpCy+YGaxSexeurijJPb9fEduYAqWUB/A+MKbQCyn1mFIqXikVn5iYWPQohSiijz4yHaR9+tzY1/fpY4ZOvvOOjJwRrq0oyf04UDvX/VrAb7nulweaA6uUUoeBaGBBfp2qWuvPtNaRWuvIKtcOOhbiJm3bBuvWwRNPgMcNjgPz8IDnnjPLFcTFOTY+IUpTUX4F4oBGSql6SikfYAiwIOdBrfVFrXVlrXVdrXVdYCPQT2sdXyIRC3EdH38Mvr5m5unNGDoU/Pzgyy8dEpYQlig0uWutM4BRwFJgNzBDa71TKfWqUqpfSQcoRFGkpJg1ZO67DypWvLlrBQfDoEFmclNSkmPiE6K0FenDq9Z6sda6sda6gdb6jexjr2itF+RzbmdptYvS9v33cOkSDBvmmOuNGAGXL8OsWY65nhClTWaoCrfwzTdmjfaYGMdc7/bbzczWb791zPWEKG2S3IXLO3sWfvjB1Mo9PR1zTaVMiefHH+HMGcdcU4jSJMlduLzZs83kpQcecOx1Bw+GzEyYM8ex1xWiNEhyFy5v7lyzw1LLlo69bkQENG4MM2Y49rpClAZJ7sKlXbwIK1bAgAGO3w9VKdN6X7kSZM6dcDWS3IVLW7zYlGT69y+Z6/fvD1lZ5nmEcCWS3IVLmzsXQkPNHqkloVUrqFkTFi4smesLUVIkuQuXlZZmRsn063fjyw0URino2xeWLpV13oVrkeQuXNbatXDlCvTuXbLP07eveZ5Vq0r2eYRwJEnuwmX98AN4e0PXroWceOWK2YmjY0eoV89s0TRmjNmZowi6dgV/fzMLVghXIclduKwffjAzScuXL+Ck5cvNGr4vvghXr5ovqFYNPvwQmjWDt94qdG1fX1/o0sWUZoRwFZLchUs6dgx27ixkt6VvvoFevcxKYOvXQ3w8fPUVLFsGBw+axdvHjoWRI82QmALceSfs31/kxr4QlpPkLlzSkiXm9rrJ/YcfzCpit98OGzbkHU5Ts6bZbumll8zavn/9a4HP16OHuV2+/ObiFqK0SHIXLumnn6B6dQjLs1U7sG+fmX3UogUsWABBQflfRCl47TV4+mlTk//qq+s+X+PGcMstUpoRrkOSu3A5WptZo1265DMrNSMDHnoIvLxMYg8MLPhiSsHEiWY5yVGj4PDh657Wo4eZDSubZwtXIMlduJy9e+HUKejcOZ8H338fNm0ym6nWrp3PCfnw8oIpU8z/H330uh2s3bub5Q5++eVGohaidElyFy5n5Upz26XLNQ+cPAn/+IeZ1TR4cPEuWrcu/Otfpmk+d26+p3TqZG7XrCnepYWwgiR34XJWrTL9oQ0aXPPA3/5mpq2+996NrSL22GMQHm7GwF+9mufhatWgSRNYvfqGwhaiVElyFy5Fa5Pc89Tbd+2C//0Pnnsun6xfRF5e5g/D4cPwxRf5ntKpE/z8s1nnXQhnJslduJTdu+H06Xzq7a+/bqaRFjKksVDdupnhk//8Z76t95gYU3ffvv3mnkaIkibJXbiUfOvte/fC9OlmSGPlyjf3BErB+PHw229m/Ps1curuUpoRzk6Su3Apq1aZQTD16uU6+P774OMDf/6zY56kSxeIjjZDJK+pv+Q8tyR34ewkuQuXkZWVT7397Fkz+ehPfzILuzuCUvD883DgQL6rhcXEmBEzhSxJI4SlJLkLl7FrF5w5c029/fPPISXFdKQ60j33mCmpEyfmeahTJ/M3Zdcuxz6lEI4kyV24jDz19qws+OQTcyA83LFP5uUFTz1lPirs3Wv3UEyMuZXSjHBmktyFy1i3DmrVMvONALOK15Ej8PjjJfOEDz9skvznn9sdrlfPjLOXyUzCmUlyFy4jz+KOX3wBlSqV3O7YoaFw991maYJce+wpZVrvq1dL3V04L0nuwiX89hscPZoruZ85A/Pnm2V9y5UruSceOdIU2BctsjscEwO//276XIVwRpLchUvYuNHc2pL77NmQnm5WgCxJ3bpB1arw7bd2h2+7zT4uIZyNJHfhEjZsMEPZW7XKPvDtt2b7vBYtSvaJvbzgvvtg4UK4dMl2uFkzs5qwJHfhrCS5C5ewYQO0bp1dgTlxwvRmDh16YwuEFdfQoabmPm+e7ZCnJ7Rta1YXFsIZSXIXTi8tzWx/aivJzJhhejKHDCmdANq3N0N0rinNtGsHCQlmmL0QzkaSu3B6CQmm4WxL7tOnm2Z848alE4BS5g/J8uWQmGg7HB1tdmWSzTuEM5LkLpzehg3mtn174OBBiI0tvVZ7jqFDzTozs2bZDrVrZ26l7i6ckSR34fQ2bDCTl2rVAr77zhws7k5LNys83OzGnas0ExpqqjWS3IUzkuQunN6GDaYEAsDMmdChg1n3pTQpZf6grF1rFpTPFh0tnarCOUlyF07t5EmzwkD79phZTFu3ltyM1ML062c6cnNNaGrXDo4dM5OshHAmktyFU7ObvJSz/G7fvtYEExFhFnRfsMB2KOcThbTehbOR5C6c2saN4O1tBsewcCE0amR2qbaCUuYPy7Jlti34WrUyk6uk7i6cjSR34dTi4qBlSyiXfgV++sm6VnuOfv0gOdnEgplU1bKlJHfhfCS5C6eVlQWbN0NkJGaMeVoa9OljbVCdO5t1BxYutB2KjjaTrDIyrAtLiGtJchdO68ABs5xLZCQmmQYHQ8eO1gZVrhz06GHiyV7vt10705jfscPa0ITITZK7cFrx8eY2snWWGaHSs6cpwFutXz+zvs2WLQBERZnDOfEK4QwkuQunFR8Pvr4QdiXWjC23ut6eo1cv07maPWqmQQMICjIlJCGchSR34bQ2bzadlV5LF5llGHv2tDoko0oVU2hfsgQADw8zmkeSu3AmktyFU8rMNAty2TpT27aFChWsDusPPXqYoTznzwPQpg1s22b2DxHCGUhyF05p3z64cgUimyWZJNq9u9Uh2eve3QznyR4S2aaNWbly506L4xIimyR34ZRsnalp60wSdbbk3rYtlC9vJjRhkjtIaUY4D0nuwinFx4O/PzTdNdck0Zz1dZ2Ftzd07WqSu9Y0bGjClOQunIUkd+GU4uNNJ6XnimXQpYtzDIG8VvfucPgwHDwonarC6UhyF04nI8Ms/hjZ6CL8+qvzlWRy9OhhbpcvB0xpJiFBOlWFc5DkLpzOnj1mxmck2YV3Z03uDRtCnTp2dffUVNi1y+K4hECSu3BCOZ2pbY7PN0vsltZeqcWllPnD89NPkJEhnarCqUhyF04nPh4CAzWNY78xyVMpq0O6vh49zAI4sbE0aiSdqsJ5SHIXTmfrVmjZ4AoeF887b0kmR9eu5nblSjw8zPruktyFM5DkLpxKVpaZ6Rnhv98c6NzZ0ngKVakStGgBq1YBf3SqyvK/wmqS3IVTOXwYLl+GiMtroWlTqFbN6pAK17kzrFsHaWlERppNmqRTVVhNkrtwKgkJ5jbi17nO32rP0bkzpKRAXJzZDhDbasBCWEaSu3AqCQng4aFpnrwJYmKsDqdoOnUyt6tW0agR+PmZfgMhrCTJXTiVhARoVOkc/qS4TnLPVXf39DT/zfkEIoRVJLkLp5KQABFeO6FJE6he3epwii5X3b1lS9Nyz96FTwhLSHIXTuPSJTh0CFqcW+069fYcOXX3+HgiIswy78eOWR2UKMuKlNyVUncppfYqpQ4opcbm8/gTSqntSqmtSqm1Sqkwx4cq3N22beY2ItWF6u05br/d3K5aRcuW5r9SdxdWKjS5K6U8gQ+BnkAYMDSf5D1Nax2utW4JvA285/BIhduzJXcSXC+5V64M4eGwahXh4WZSrSR3YaWitNzbAge01r9qrdOA6cDduU/QWl/KdTcAkGqjKLaEBKjgfZlajfyhRg2rwym+7Lp7oE8ajRpJp6qwVlGSe00gd/XwePYxO0qpp5VSBzEt92cdE54oSxISNBF6K6qzi7Xac3TqZJaz3LLF1qkqhFWKktzzW7UpT8tca/2h1roB8Ffg5XwvpNRjSql4pVR8YmJi8SIVbi0zE7Zv00RkbIaOHa0O58Z06GBu160jIsIsRX/xorUhibKrKMn9OFA71/1awG8FnD8d6J/fA1rrz7TWkVrryCpVqhQ9SuH2Dh6E5BQPU2/PSZKupnp1qF8f1q2zdarm9CMIUdqKktzjgEZKqXpKKR9gCLAg9wlKqUa57vYG9jsuRFEW2JYdqHjcJEhX1aGDSe4R5sOtlGaEVQpN7lrrDGAUsBTYDczQWu9USr2qlOqXfdoopdROpdRW4M/AsBKLWLilhATwJIOw2ys59/rthenQAU6donryQapUkU5VYR2vopyktV4MLL7m2Cu5/v+cg+MSZUzCpqs05QC+ndpaHcrNyS4pqfXraNmyobTchWVkhqpwCglbsmjBNtett+cIC4OQEFvdfccO2TBbWEOSu7DcuXNw7Kw/EV67zFZGrszDA267zTZiJjUV9u61OihRFklyF5bbvt3cRjRLAx8fa4NxhA4dYNcuWtYz4yCl7i6sIMldWC4hNhWAiC4VLY7EQbJLS00S11KunIyYEdaQ5C4sl7DyLFU4TbUeLawOxTGiosDbG6+NawkPl+QurCHJXVguIcEsFqZua291KI7h7w+tW9s6VWVtd2EFSe7CUhkZsOP3ykRUOgEVKlgdjuN06ACxsUTcms6ZM/BbQXO6hSgBktyFpfbtySI1y4eI8CyrQ3GsDh0gNZWWvmaojHSqitImyV1YKmHRcQAiulW2OBIHy+5UbXHmJ0Dq7qL0SXIXlkr46QzepNF0YHOrQ3Gs0FBo2JCgzStp0ECSuyh9ktyFpRK2e9LMaz8+TepZHYrjZS8iFhGhJbmLUifJXVgq4XQ1Imqcce3Fwq6nQwdITCSi5lkOHICkJKsDEmWJJHdhmcQdpziZGUpECzcdJ9jeDO2M8NiO1n/MxBWiNEhyF5bZNucA4IadqTmaNYPy5Yk4azpVZeMOUZokuQvLJKy5AEDEPQ0sjqSEeHpCVBR1dv1AUJAMhxSlS5K7sEzCTm+qeydS5RY/q0MpOdHRqG0JtGieKcldlCpJ7sIamZmmM7X6aasjKVnR0ZCRQUS1U2zbJssQiNIjyV1YIm3bHnZlNSUi3M2zXbt2AER47ODyZTh82NpwRNkhyV1YYs/8vaTjQ4vObrLM7/VUrQr16hFxbiUgdXdReiS5C0skrM7uTL2zmsWRlILoaG7dMxulJLmL0iPJXVgiYYcn5TzSaNKsDPwIRkcT8Nt+GtVLl+GQotSUgd8s4XSuXCHhTE1uDT2Dl5fVwZSC6GgAIkJPSctdlBpJ7qLU6bh4EmhBRPNMq0MpHRER4ONDC88dHDwIly9bHZAoCyS5i1J3asUOEqlKRIwbbc5RkHLloHVrIs6tAmQZAlE6JLmLUpfw01kAIjoEWhxJKYqOJuLgHECWIRClQ5K7KF1ak7DDEzDVijKjXTtqp+4npHyG1N1FqZDkLkrX8eMkXK5H7QqX3WrL1EJFR6OAFqGnJbmLUiHJXZSuTZtIIIKIWzOsjqR01akDoaFEeO1g2zbIcrMtY4XzkeQuStXVtfHsoSktOgRZHUrpUsrU3c+vJikJDh2yOiDh7iS5i1K1a3UimXgR0drT6lBKX7t2RJxaCshMVVHyJLmL0pOeTsJOM2upTHWm5oiO5lZ24uGhJbmLEifJXZSeHTtISG+Gn08GDRtaHYwFIiPx80ijccUzktxFiZPkLkrPxo0kEEF4sww8y2BVhvLloXlzIrx2ylh3UeIkuYtSozduIkG1JCKqnNWhWKddOyIurOHQIbh0yepghDuT5C5KzYl1hzmvKxDRUlkdinWio2lxdRMgM1VFyZLkLkrHhQskHAwAymhnao7oaCIwBXdJ7qIkSXIXpSM2lgRMVg8PtzgWKzVtSs3yl6lY7op0qooSJcldlI7smal1b8kiONjqYCzk4YFq15YW3rsluYsSpbRF27FHRkbq+Ph4S57bWVy6dInTp0+Tnp5udSgl7/RpfkupiJefF1WrWh2MxS5c4NxFD66oIGrXNpNXhfD29qZq1aoEBRU8e1sptVlrHVnY9crCPjhO6dKlS5w6dYqaNWvi5+eHcuffcK3JTE0nKaA51asrata0OiCLXbjAmQPnOUw96tcHX1+rAxJW01qTkpLCiRMnAApN8EUhZRmLnD59mpo1a+Lv7+/eiR0gNZWrGd6Awt/f6mCcQEAAfqQAkJxscSzCKSil8Pf3p2bNmpw+fdoh15TkbpH09HT8/PysDqN0JCWRjHmtZeUlF8jbGz+fLECTkmJ1MMKZ+Pn5OaxMK8ndQm7fYs+RlESK8sfDQ1OuDM9fys0j0B9flSotd2HHkTlBkrsoeUlJJHsE4u+vpPMwR0AA/jqJlBRrBjQI9yfJXdyw8ePHo5Sy/fP39yc8PJzPPvvsj5OystDJyaRk+ZZYSSYnjkaNGuX7eMOGDVFKMX78+GJd9/Dhwyil+P77723H6tatywsvvHAz4RrZdfe0NEVGMfctyS+u0tC5c2cGDhxou79s2TImTpyY57yHH36YyMhCB3OIEiajZcRNCQ4OZsmSJQAkJSWxcOFCHn/8cQIDA7n//vshOZk07U0mHiXamerr68uhQ4eIj4+3SyxxcXEcOXIEX2cbkuLvjz8nAUhJMWuKObuPPvoIb29v2/1ly5Yxa9YsRo8ebWFU4nokuYub4uXlRXR0tO1+t27dWL9+PfPmzTPJPSmJFExWL8nO1ICAAFq3bs306dPtkvv06dPp2rUrmzdvLrknvwFX09Lw89OQYkbMOHNyT0lJwc/Pj7CwMKtDEcUgZRnhcOXLl/+jxz+73n7ixCEeeKA/QUFBlC9fnr59+3LgwAHb18ycORMPDw9WrFhhO3b48GGCgoJ4+eWXi/S8Q4YMYcaMGeRMzNNaM2PGDIYMGZLn3GtLDACrVq1CKcWOHTuK9XrXrl1LTEwM/v7+VKpUiZEjR3L58mXb41OmTEEpRWxsLJ07d8bPz4933nkH7/K+eJFuV3c/efIkw4cPp379+vj5+dG4cWNefvll0tLSCowhNTWVJ598kpCQECpVqsSLL77IxIkT83TQHTp0iP79r/8+gOnUe++99xg9ejRVqlQhPHu9iNzfs/Hjx/Puu+9y5MgRW1nu4YcftrvO8uXLadGiBQEBAXTs2JGdO3fmeZ7333+fMWPGUKlSJSpXrsyECRMA+N///kf9+vUJCQlh+PDhXL16tQjvhMhNWu7ipmVkF42Tk5NZsGABq1ev5r///a95MCmJ85nBPP30bQQEePP555/j5eXFuHHjiImJYfv27VSsWJFBgwYxZ84chg8fzvbt2ylfvjzDhw+nXr16vPLKK0WK45577uHJJ59k7dq13H777fz8888kJiYyYMAAXnzxxRJ57evWraNbt27079+fWbNmcfbsWcaOHcv58+eZNWuW3blDhw7lySefZNy4cYSEhKCy6+4pSQGAWeD+zJkzVKxYkffee48KFSqwb98+xo8fT2JiIp9++ul14/jLX/7ClClTePPNN2nWrBmTJ09m+vTpduekpqbSrVs3vL2v/z7keOedd+jUqRNff/01Wfns5v3oo4+yf/9+fvrpJ+bOnQtAlSpVbI8fPXqUF198kZdeegk/Pz9eeOEF7rvvPnbs2GH3B+fdd9+ld+/efPvtt3z//fe8+OKLnD59mri4OCZNmsTRo0d5/vnnady4MWPHji36GyNM68aKf23atNFl2a5du6wO4aaNGzdOA3n+Pfvss+aEtDSt4+L0y3/7j/b09NQHDx60fe2xY8e0t7e3fvPNN23Hzp49q6tXr66HDx+u//3vf2tvb2+9devWIsVRqVIlrbXW/fr100899ZTWWusnn3xS33333VprrStVqqTHjRtn+5qYmBh977332l1n5cqVGtDbt2/XWmt96NAhDeiFCxfazqlTp44eM2aM7X7Hjh11586d7a6zYsUKu+tMnjxZA3rixIn2gaek6KNxJ/Xm+CydlZX/a0tPT9dTp07V5cqV06mpqfnGdebMGe3r66vffvtt29dlZWXpsLAwbX7FjY8//rhI7wOgW7ZsmSeWa79nY8aM0XXq1Mlz3rBhw7Snp6fet2+f7djcuXM1oHfv3m33PLm/d5mZmbpatWo6JCREX7x40XZ80KBBum3btvl/g9xQYbkBiNdFyLHScncmo0fD1q3WPHfLlpDPyIfCBAcH8+OPPwKmZbh582ZeeeUVKlasyLjnniMTDxK2byY8vDX169e3fV2tWrXo0KEDa9eutR2rWLEin3/+OX369MHHx4dx48YRkWt94MzMTFvJBcDT0zNP2WHIkCGMHj2a9957j1mzZjFp0qRiv6aiSk5OZsOGDXzwwQe2Ty8AHTt2xNvbm82bN9O8eXPb8d69e9tfoFw5/DxSycpSXL1q+iS01vz73//ms88+49ChQ3bliKNHj9Iwn/0Jt2/fztWrV+nXr5/tmFKKvn37smvXLtux2NhYWrcu/H3IN9Ziqlu3rt3opZx6/fHjx2natKnteLdu3Wz/9/DwoF69evj7+9tNv2/YsCHr16+/qXjKIqm5i5vi5eVFZGQkkZGRdOjQgWeffZa///3vvPnmm5w7fpwU/Dlz5iShoaF5vjY0NJRz587ZHevatSuhoaFkZWUxcuRIu8dySgo5/1avXp3nmv369ePKlSu89NJLJCUl0bdvX8e+4FzOnz9PZmYmTz31lF1c5cqVIz09nWPHjtmdn+d7oBT+/uaPU85M1YkTJzJmzBgGDBjA/PnziY2N5cMPPwS4bt35999/B+zLIvndP3my6O9DfucVR0hIiN19Hx8fIO9ryO+8/I5Jzb34pOXuTG6g5eyMwsLCSEtL4+CePdRpdBuVK1fn99935jnv1KlTdnVegLFjx5KZmUm1atUYPXo006ZNsz326aef2nVUNmnSJM81AwIC6NOnD++//z6DBg0iICAg3xh9fX3zdFJem+AKExISYhs/36tXrzyP16hRw+5+frMPfct7oa5kkZIEVPRg5syZDBo0iDfeeMN2Tu7Wd36qVasGQGJiot33MzEx0e686tWr5+nUhPzfhzIze9qNSXIXDpcz2qR2cDApnuUJD2/H4sVfcejQIerVqwfAiRMnWL9+vd3EotWrV/PBBx8wY8YMgoKCuPPOO7n33nu59957gfyTeX6efPJJUlNTeeKJJ657Tq1atVizZo3dseXLlxfnZRIQEEB0dDR79+4tcqfvtTwCA/DlKslXvAEPUlJSKHfNGg1Tp04t8Brh4eH4+voyf/58/vKXvwCmvLNw4UK789q1a8dXXxX+PhSVtKidmyR3cVMyMjLYuHEjAGlpaWzevJnXX3+du/v0oVqFCuzWvgwe/DDffvsWPXv25NVXX8XT05Px48dTuXJlHn/8cQCuXLnCI488wuDBg23D7R5//HGefPJJOnXqlKfEUJDOnTvTuXPnAs8ZMGAAX375Jc8//zy9e/dm5cqVLF26tNiv/+2336Zbt254eHgwcOBAypcvz9GjR1m0aBFvvPEGjRs3LvgCAQH4cZErV01C7969O5MmTaJdu3Y0aNCAqVOn5hmqeK2c4Zfjxo3D29vbNlrm0qVLdi3whx9+mLfeKvh9KI6mTZty6tQppkyZQvPmzalcuTJ169Yt9nVEyZCau7gpFy9epH379rRv35477riDjz76iCeeeIKvP/gADaSkexESUo4ff/yRpk2bMmLECIYNG0adOnVYtWqVrRwwZswYUlJS+M9//mO79oQJEwgMDCywBX6jevfuzZtvvsmsWbMYMGAAR44cyXcqfWE6duzImjVrSExM5MEHH6Rv3768/fbb1K5du2h1ay8v/L3SSMv0JCMDXnnlFYYOHcrLL7/M0KFD8fHxKVKn8Ntvv83DDz/M+PHjGTp0KKGhoYwYMcKuY7JcucLfh+K47777ePjhh/nLX/5CVFTUDbX+RcmRnZgssnv3bpo1a2aepBfHAAAgAElEQVR1GCXnyBGunr3CjqxbqVMHitHwLnMu7v2d/Zer0aSxpnyQ42rdd9xxB+np6fl2PAvnVVhukJ2YhLWSkkgpVwFSkA06CuEf5AWXIflKJuWDbuxXcuXKlWzatInWrVuTnp7Od999x4oVK5g5c6aDoxWuQpK7cLzMTEhOJjmwFiAbdBTGO8gPrxPppFzO4kZ/JQMDA5k3bx7//Oc/uXr1Ko0aNWLKlCl5llgQZYckd+F42TtQJGs/fH3BQ3p2Cubnhz9XSL564ytXRkVF2Tq2hQDpUBUlISkJgJQ0LynJFIWHB37e6aSke2FRF5hwQ5LcheNduUK6jz9p6bIhdlH5+4HGg6speRfpEuJGSHIXjpfTmYp0phaVX3ZHaspFx2yOLESRkrtS6i6l1F6l1AGlVJ51N5VSf1ZK7VJKbVNKrVBK1XF8qMIlpKVBejrJnoGAdKYWlW+IL4oski8Xc889Ia6j0OSulPIEPgR6AmHAUKXUtVuybAEitdYtgFnA244OVLiIK1cASM7yxdsbcu3KJgrgUc4HX5VKSoqs6SIcoygt97bAAa31r1rrNGA6cHfuE7TWK7XWydl3NwK1HBumcBlJSaCUdKYWl1L4e6eTnCF/DYVjFCW51wRyr116PPvY9YwAfriZoIQLS0oiyz+QlKvSmVpcfn6Qrr1JvyqlGXHzipLc8/ucmO+ALaXUn4BI4J3rPP6YUipeKRV/7XKkwg1kZWV3pgYDrtuZeuXKFZRSTJkypVSf1798dqfq+dRSfd5rDRw4sNCF1xypbt26vPDCC6X2fIBtwbQcOdsZXrhwwe68nP1vr2SXG11JUZL7caB2rvu1gN+uPUkpdQfwEtBPa53vT6fW+jOtdaTWOrI4q/wJF5GSAlqT7FkekM7U4vILMStDpkinaol79NFH7VYB3bdvH//4xz/yJHdXVpQZqnFAI6VUPeAEMAS4P/cJSqlWwKfAXVrr0w6PUriGnMlLWb54eMA1y5LftMzMTDIzM227+rgbb19PvFU6ydfpVE1JScFP/mLelPT0dDw8PKhVqxa1arl312ChLXetdQYwClgK7AZmaK13KqVeVUrlbNr4DhAIzFRKbVVKLSixiIXTePjhh4mMjGT58uW0aNGCgLp16ThyJAk7duPvDzlLiSuleO+993juueeoWLEiISEhPPPMM3l2Qrre9efNm8ett96Kr68vmzZt4uTJkwwfPpz69evj5+dH48aNefnll+2ud/jwYZRSzJgxg8cff5zg4GBq1arFuHHjyMqynyg0e/ZsGjdujJ+fH506dWLPnj15YsnMzGT8+PHccsstlCtXjltvvdVul6jc8S5atIiwsDD8/f3p3bs3586d48CBA3Tp0oWAgAAiIyPZtm1bvq/Zz8vMVF21ciVKKZYuXUq/fv0IDAxk1KhRALz77rtERUURHBxMaGgoffv2zbPme+fOnRk4cCDTpk2jYcOGBAUF0bNnT44fP2533rFjx+jVqxd+fn7UrVuXL774Ik9MOSWMLVu2EB0djb+/P61ateLnn38u4N0zxo4dS3h4OIGBgdSqVYsHHnjAti1gQf7zn/9Qu3ZtAgIC6N+/PytWrEApxapVq2znJCcn8+yzz1KtWjV8fX2Jiopi2bJl+X4fPvvsMxo0aICvry+//fabXVlm1apVtu0Y69Wrh1Iqz7r0hw4donv37gQEBNC0aVPmzJmT7/NMnjyZevXqERgYyIMPPkhqaiqxsbG0bduWwMBAOnfuzNGjRwt9/Q5RlF20S+JfmzZtircluJspbIdzVzBs2DBdpUoVHRERoadPn67nT5qkG9Wtq+vXD9OHD2fZzgN0jRo19D333KMXL16s33nnHe3j46NfeOGFQq9fqVIl3ahRI/3111/rZcuW6WPHjult27bpMWPG6Llz5+pVq1bpzz77TNeoUUM/9thjtq89dOiQBnSdOnX0n//8Z71s2TL917/+VQP6u+++s523efNm7enpqQcOHGiLrV69ehrQkydPtp33t7/9TXt5eenXXntNL1myRI8cOVIDetq0aXm+H61bt9azZ8/WX3/9tQ4JCdH33nuvbtOmjf7kk0/04sWLdUREhG7WrJnOyvrje5Tj2N4kHR+XqVcsWaoBXbNmTf3yyy/rFStW6HXr1mmttR49erSeMmWKXrlypZ4/f77u2bOnrlq1qr5w4YLtOjExMbpWrVq6ffv2et68efrbb7/VVatW1T179rSdk5WVpVu1aqVr166tp06dqmfPnq2bN2+ua9SooWNiYmznjRs3Tvv5+enw8HD93//+Vy9evFi3a9dOV6pUSSclJRX4Hj7yyCN62rRpetWqVXrmzJk6OjpaN2vWTGdkZNjOqVOnjh4zZozt/pw5czSgn3rqKb106VL997//XdepU0cDeuXKlbbz7r//fh0YGKgnTZqkFy9erAcMGKC9vLz0zz//bPd9qFatmm7ZsqWeOXOmXrRokb548aIeN26crlSpktZa64sXL+oJEyZoQM+ZM0dv2LBB//LLL1prrSdPnqwB3bx5cz1p0iS9dOlS3adPH+3t7a2PHTtm9zw1a9bUMTExeuHChfrDDz/UPj4+euTIkbpFixb6m2++0XPnztW1a9fWd955Z4Hfs8JyAxCvi5BjJblbJL838LnntI6Jsebfc88V/zUMGzZMe3p66n379mmdnq51XJz+7nPzy7B+/W7beYBu0qSJzszMtB17/fXXtZ+fnz579myB1wf0li1bCowjPT1dT506VZcrV06npqZqrf9I7g8++KDduREREXrw4MG2+4MGDcqTaF9//XW75H727Fnt7++vx48fb3etnj176saNG+f5fhw4cMB27MUXX9SA/t///mc7tmjRIg3k+zNw5rerOi5O6yWzFmpAjx49usDXnpGRoZOTk3VgYKDdc8TExOigoCB97tw527H3339fAzo5Odkujo0bN9rOOXz4sPb09MyT3AG9YsUK27EtW7ZoQP/www8FxndtrMePH9eAXr16te34tck9MjJS9+rVy+5rn3zySbvkvmvXLq2U0lOmTLGdk5mZqW+99Vbdo0cPu++Dr6+vPnnypN31cid3rbVeuNB8vw8dOmR3Xk5y//LLL23Hzpw5oz09PfXHH39s9zzBwcF2f2AHDRqU57V++OGHGijwj6KjkrssPyBuSt26dWnUqJGt3l63SWsAzp61//h/991345Frech77rmHlJQU236r11OzZk1atmxpd0xrzcSJEwkLC8PPzw9vb28eeOABUlNT83zk7dGjh939sLAwu9JEbGws/fr1s9uO7p577rH7mh07dpCcnMygQYPsjg8ePJh9+/Zx+vQf3Ux169alQYMGtvsNGzYEoGvXrnmOnThxIs/r9Q82/QlXkzMBs2PUtTZu3Ej37t2pVKkSXl5e+Pv7c+XKFfbt22d3XlRUFBUqVLB77bmfNzY2ltDQUNq1a2c7p06dOrRp0ybPc3p7e9uNoMm51rVlnmv98MMP3HbbbQQHB+Pl5WWrc18ba47MzEy2bt1Kv3797I5fez8uLg6ttd174uHhwaBBg1i7dq3duW3atLFtIn6jcv8cVapUiapVq+Z57ZGRkQQHB9vuN2zYEB8fHzp27Gh3DOC33/KMSXE4WfLXidzALm+WCwkJMf/JTu5Z2csOZGXZb5xctWrVfO+fPHmywOvnt1XdxIkTeeGFFxg7diwxMTFUqFCBuLg4nn766TwbNtviy3btps6///77dWPLkRPjtbHk3D9//rzta/J7vmuP5xzLb3Ppcr4KRRZpqfbPkePo0aP06NGDtm3b8umnn1KjRg18fHzo3bt3kV577ufN77XnvP7Lly/bHQsKCrL741zQa8gRFxdHv379GDBgAGPHjqVq1aoopYiOjr7u1yUmJpKRkZFnz9xr7588eZLAwED8rxlvGxoaSnJyMqmpqbaNxou03WEhCvs5ut455cuXL/b3zVEkuQvHuHIF/PxIuWx+kK9dwz136zb3/erVqxd42dwt6hwzZ85k0KBBvPHGG7Zju3btupGoqVat2nVjy5ET4+nTp6lUqZLt+KlTpwBuaP/R6/HwAD+vDK5megJ5X/+SJUtITk5m/vz5BAQEAGaT8nPnzhX7ufJ77WBepyNG5cydO5cqVarw3Xff2V7HkSNHCvyaKlWq4OXlxbXzYK69X716da5cuUJycrJdgj916hT+/v62xA75/wyVBVKWETdPa9NyDwggJSX/U+bPn283SmXOnDn4+fnRvHnzYj9dSkqK3S8vwNSpU4t9HTCliwULFpgOqFyx5da8eXP8/f3zbFk3Y8YMGjdunKdVebP8/DSp5D+ONCUlBQ8PD7y8/miXzZgxg4yM4o+Nj4qK4tSpU2zatMl27OjRo/zyyy/FD/o6sXp7e9sl18LeJ09PT1q2bMn8+fPtji9YYD8ALyoqCqUUs2bNsh3TWjNr1iy7MkhRlWaLurRIy13cvNRUyMwk3a88GRmX8z3l8uXLDBo0iJEjR7Jz505effVVRo0aZWv1HjlyhAYNGvDf//6Xhx56qMCn6969O5MmTaJdu3Y0aNCAqVOn5hkKWFR//etfadeuHffddx8jRoxgx44dfPnll3bnVKxYkdGjR/P666/j5eVFZGQkc+bMYfHixXz77bc39LzX6tatGwArVqzAv7wnmeS/xkzXrl3JzMzkkUceYcSIEezcuZMJEybkKQkURa9evYiIiGDQoEG89dZb+Pr68sorr+RbqrkR3bt3Z+LEiYwePZq+ffuyfv16vvnmm0K/7m9/+xv33HMPo0aNol+/fqxbt45FixYB2EoczZo1Y+jQoYwaNYpLly7RsGFDPv/8c/bs2cPHH39c7FibNGkCwKeffsqQIUPw9/cnPDy82NdxJtJyFzcvZyVIj8DrnjJmzBiqV6/O0KFDefXVV3n00Ud58803bY9rrcnMzMwzBj0/r7zyCkOHDuXll19m6NCh+Pj4MGnSpBsKPTIykunTp7Nlyxb69+/PvHnz+O677/Kc9+qrr/J///d/fPzxx/Tp04c1a9bwzTffMGTIkBt63mvlTNCCP5YhyE94eDiTJ09m06ZN9OnTh2nTpjFz5ky7jryiUkqxYMECwsLCGD58OKNHj2bUqFG0b9/+hl9Hbr169eKtt95i9uzZ9OvXj9WrV/P9998X+nUDBgxg0qRJzJs3j/79+xMXF8eECRMAU/vP8fnnnzNs2DBee+017r77bo4cOcL3339/Qy33OnXqMGHCBObMmUOHDh1s495dmcr9cbQ0RUZG6vj4eEue2xns3r2bZs2aWR2GYxw5AufO8Xu1lhw/oWjZEnJVDVBK8cEHH9gm4YiCZWTA1q1Q0/N3qre6uVEe7uL111/njTfe4Ny5c24/S7ew3KCU2qy1jizsOlKWETcvu96enKLw8bFP7KL4vLygnFcGyRk+kJ5e5hbFT0xM5J///CddunTB39+fn3/+mbfeeosRI0a4fWJ3JPk1FDcnMxOSk6F6dZLPu+5KkM7G30+TfDnA/OG8gXq6K/Px8WHPnj189dVXXLx4kerVq/Pcc8/x2muvWR2aS5HkLm5O9vj2TP/yXD0JuebM2FhV+nNl/uU9OX/Zm4xLZ/AqY8k9ODiYxYsXWx2Gy5MOVXFzsjtTUzzMmOvsodfiJgUEml/N5MuZFkciXJUkdwu5RYs2e/JS0lUz6UbKMo6R831MTvEwm6CIMsGROUGSu0W8vb1Jud6MH1eRM3kpMJDkZNPv56ZLrZc6Ly/w8cokCX/TpyHKhJyJX44gyd0iVatW5cSJEyQnJ7tuCz4lxXSoBgSQnCytdkfzD1Ak428rfQn3pbUmOTmZEydOOGwSmXSoWiRnMsZvv/1Genq6xdHcoMuX4dw5srzLcexkIsHBZoy2cIyLF+HCBfBIOoPH+fNWhyNKmLe3N6GhoXYTtW6GJHcLBQUFOeyNtMRDD8HSpayb/Ts9eyrmzwcHTW4UwA8/QK9esLpCfzqdnfvH1lZCFIGUZcSNW78eOnRg8y8m6eSzDLi4Ca3N0vhsPl8Pfv3V2mCEy5HkLm7MqVNw8KBJ7pshNBRq1LA6KPcSGgo1q6bzC63NH1IhikGSu7gxOcnmttvYvBkiI6VqUBJat/XiF49IWLfO6lCEi5HkLm7MunVQrhxJTVqze7eUZEpK6zaKPVmNSVq7xepQhIuR5C5uzPr1EBlJwp5yZGVJci8pbdpAFp5s2eljhs4IUUSS3EXxXb0Kmzfb6u0gyb2kREWZ2zgiIdeOSUIURpK7KL7NmyEtzVZvl87UklOtGtSqmUUcbaXuLopFkrsovpwkk53c27SRztSS1LadB3E+HWTEjCgWSe6i+Navh8aNSQ6owq5dUpIpaVFRcCDtFs5t3CdTgEWRSXIXxaO1Se7ZrfasLGjb1uqg3FtO3T0+qSkkJFgbjHAZktxF8ezZA4mJcPvtxMaaQ5LcS1bOJ6M4omDNGmuDES5DkrsontWrzW1MDLGxULcuOGgRO3EdISHQpAnE+neW5C6KTJK7KJ41a8zQmPr12bRJWu2lJSoKYnUUevUa2bxDFIkkd1F0WpuWe0wMp04rjhyR5F5aoqLg95QQTpz3g127rA5HuABJ7qLofv0VfvsNOnUiLs4catfO2pDKipw/onFE/VEaE6IAktxF0eWqt2/aBJ6e0KqVtSGVFS1bgpeXJq58N6m7iyKR5C6KbvVqqFIFmjYlNhaaN4eAAKuDKht8fSE8XBEb0Nkkd1fdmlGUGknuoujWrIFOncjSithYKcmUtqgoiL/YiKzfT8H+/VaHI5ycJHdRNEePwuHD0KkTBw6YBQqlM7V0tW0LF1PKsZ9GUpoRhZLkLoomJ5lkj28HabmXtpzv98agO6VTVRRKkrsomtWrzWya8HA2bTK19mbNrA6qbAkLg6Ag2FC5r7TcRaEkuYuiWb0abr8dPDyIjTXb6nl6Wh1U2eLhYVrvG662+qNMJsR1SHIXhTt50nTgxcSQmgpbt0pJxiq33QbbT1biEuWl9S4KJMldFC4niXTqxC+/mH06JLlbo3170FoRW/4OqbuLAklyF4VbswYCA6FVK9t+EbfdZm1IZVXOH9UNNQdKy10USJK7KNzq1dChA3h5sW4d1K9vtn8TpS8kBG69FTbQHg4cMMtBCJEPSe6iYImJsHMnxMSgtdlhr0MHq4Mq29q3hw0napOFgpUrrQ5HOClJ7qJgP/1kbrt25eBBOH1akrvV2reHC5e92BvcDlassDoc4aQkuYuCrVhhBle3aWOrt0tyt1ZOf8eGhg/Cjz/KOjMiX5LcRcFWrIAuXWz19uBgM5lGWKdxY6hQAdaX6wLHjpnauxDXkOQuru/QIbOGe7dugKm3t29vJtMI63h4ZNfdT9c3B6Q0I/Ihv6bi+nKSRrdunD9v+lWlJOMc2reHXQfKcaFGmCR3kS9J7uL6VqyA6tWhWTM2bDCHZHy7c2jf3txuaDbcjJiRfVXFNSS5i/xpbUbKdO0KSrF+vVlLRmamOofoaPDygp/9e8DZs5CQYHVIwslIchf527HDjHvMVW9v2VJ2XnIWAQFm8bbVJxubA1KaEdeQ5C7y9+OP5rZbN9LTYdMmqbc7m06dIC6hHMmNIiS5izwkuYv8rVgBDRvCLbewZQukpEi93dnExEB6OmwMG27WmUlNtTok4UQkuYu8UlNNJ92ddwKwapU5HBNjXUgirw4dzLDINX53QnIyrF1rdUjCiUhyF3mtXWuSxV13ASa5N2smi4U5m+Bg0w+y+kQD8PaGpUutDkk4EUnuIq8lS0yy6NyZjAz4+Wfo3NnqoER+OnWCjXFepN7WxbxvQmST5C7yWrrUbKkXGMjmzXDliiR3ZxUTA1evQnzYQ7B9uywBLGwkuQt7J06YJJGrJANSb3dWHTua29WeXc1/pDQjsklyF/ZyksM19fbQUOtCEtdXubLZvGPN/mqmU0RKMyKbJHdhb8kSqFEDmjcnPd30rUpJxrnFxMC6dYqM7j1h+XLIzLQ6JOEEJLmLP2RkmMlLd90FSvHLL6be3qWL1YGJgnTqZN6nXxoPgfPnIS7O6pCEE5DkLv4QF2eSQ3ZJJmcHN6m3O7ec92dlegdQSkozAihicldK3aWU2quUOqCUGpvP452UUr8opTKUUgMdH6YoFYsXm1kxd9wBmHp7WBhUrWptWKJg1apBeDgsXxdgVnZbtMjqkIQTKDS5K6U8gQ+BnkAYMFQpde1ePEeBh4Fpjg5QlKKFC83wiwoVpN7uYu64w7xfKXcNgPh4GRIpitRybwsc0Fr/qrVOA6YDd+c+QWt9WGu9DZBFpV3VkSNm2di+fQGTH5KSJLm7iu7dzaoRP1cbZA5I673MK0pyrwkcy3X/ePYx4U6+/97c9usHmBGRSpnl3IXz69QJfHxg+f66UKeO+RQmyrSiJHeVz7Eb2m5dKfWYUipeKRWfmJh4I5cQJWXBArPzcmOzPvjSpRAVBZUqWRyXKJKAALNq5/Iflfn09eOPZilPUWYVJbkfB2rnul8LuKGCntb6M611pNY6skqVKjdyCVESLl0yQ2OySzLnz0NsrG1RSOEiunc3lbVTtw80iV3WeC/TipLc44BGSql6SikfYAiwoGTDEqVq2TKzMHh2SebHH82WnJLcXUv37uZ2xdUOUL68+TQmyqxCk7vWOgMYBSwFdgMztNY7lVKvKqX6ASilopRSx4FBwKdKqZ0lGbRwsIULoUIF224cS5ea5WRlv1TX0rq1eRt/XOVl/jJ//71snF2GeRXlJK31YmDxNcdeyfX/OEy5RriajAwzsqJ3b/DyQmuT3Lt1MxswC9fh6Wnet+XLQb/eFzVrFvzyi9lsVZQ5MkO1rFu3Ds6etdXbd++G48dtk1SFi+ne3bx/e5v0Ndl+7lyrQxIWkeRe1s2aBb6+0KsX8MeikFJvd005dfelmyqYRYFmzgR9Q4PbhIuT5F6WZWXB7NnQsycEBgImuTdtCrfcYnFs4obUq2eWaF64EBg4EPbvhx07rA5LWECSe1m2YQOcPGmSAGb03OrV0mp3dX37mvfxUrcBZq2gWbOsDklYQJJ7WTZrlpnW2KcPAGvWmC3bJLm7tr59TT/50i1VzdRVSe5lkiT3siory/zS33knBAUBMH++meko68m4tvbtoWLFXKWZXbvMP1GmSHIvq+LizLCK7JJMVpZJ7nfeCX5+Fscmboqnp+kfX7wYMvsNMAdnz7Y2KFHqJLmXVTNngre3bVbq5s1mldj+/S2OSzhE375mhOuGIzWgQwdJ7mWQJPeyKKckc8cdEBICwLx5psXXu7fFsQmHuPNOMwlt4UJg0CCz6Mzu3VaHJUqRJPeyaN06s377/ffbDs2fb/reKla0MC7hMMHBZvu9hQuBwYPNqJmpU60OS5QiSe5l0TffgL+/rQZz4ADs3CklGXfTt69prB9MqmY+pU2dKhOayhBJ7mVNairMmAEDBtgmLs2fbx66++4Cvk64nOwRrmZxyD/9CQ4fhvXrrQxJlCJJ7mXN4sVw4YL5Zc82bx60bGk28BHuo0EDiIgwfecMGGA+rX3zjdVhiVIiyb2s+eYbCA01H9OB06dNCV5KMu5p8GAzEfnouUDzJn/3HaSlWR2WKAWS3MuS8+fNGt9Dh9rW850/35RhpSTjngYPNrczZgAPPGB+Bn74wdKYROmQ5F6WzJxpWm25SjLTpkGTJubju3A/9eub5dy/+w6zZGSVKvD111aHJUqBJPeyZPJkCAszW/ZgJqiuXm0adCq/bdCFW7jvPoiPh4NHvc2bvWCBqccJtybJvazYtg02boSRI22ZfPp0U5IZOtTi2ESJuu8+cztjBub9T0+H//3P0phEyZPkXlZ8/jmUKwcPPmg7NG0atG0LDRtaGJcocXXqQHR0dmkmLMwsR/D55zLm3c1Jci8LkpNNnXXgQKhUCTCTW7ZssZukKtzY4MFmBYK9ezGt9/37TU1OuC1J7mXBzJlw8SI89pjt0LRpZkZ6zmgK4d4GDTLVuOnTs+8EB5vWu3BbktzLgs8+M0Nibr8dMJ/Gp02Dbt2gWjWLYxOlomZN6NoVpkyBLF9/M2Jq9myzdKRwS5Lc3d3OnWbKea6O1NhY+PVXKcmUNY8+alYg+OknzM9Daqp0rLoxSe7ubtIk8PWFYcNsh7780mzIMWCAhXGJUte/P1SoAF98gZnY0KED/Oc/kJlpdWiiBEhyd2dnzsBXX5kRMpUrA3DpkinJDB1qyq6i7PD1NdWYuXOzqzHPPw+HDv2xcpxwK5Lc3dknn5gdr0ePth2aNg2SkuDxxy2MS1hmxAgzSXnqVExTvm5deP99q8MSJUCSu7tKTYUPPzRb8oSFAaYj9ZNPzAqQUVEWxycsERFhliP44gvQHp7w7LOwdq2ZwirciiR3d/Xdd/D77/DnP9sOxcaasc5PPCHLDZRlI0bA9u3Z+XzECChfHiZOtDos4WCS3N2R1uaj9q23msWisn3yidmfQ0bJlG1Dh5oO9U8/BYKCYPhw0xg4ftzq0IQDSXJ3Rz/8AFu3mlZ7dhP9wgXz+3v//aahJsqu4GB46CGztP/p0/zRJ/P225bGJRxLkru70RrGjzcdZbnWkZkyBVJSpCNVGKNHm26Zjz7C/KwMG2Ymu508aXVowkEkububJUsgLg7+9jfw9gbMIoDvvQcdO9pW+xVlXNOmZo/Vjz4yf/T5v/+DjAxpvbsRSe7uRGv4xz/gllvsJi19+y0cO2Z+f4XIMWYMJCZmb6vaoIEZBP/JJ6YjXrg8Se7uZNky2LTJtNp9fADIyoK33oLwcOjZ0+L4hFOJiYFWrUzfe1YW8NJLZhD8u+9aHZpwAEnu7iIrC155BWrXhkcesR3+/nvYtQvGjpXhj8KeUqb1vnu3qebRqJHpcf/wQxk54wYkubuL774zA9lffdXWatca/vUv01+WsxCJzukAAA4mSURBVBuPELnddx/UqgVvvJG9d8err5q1Zl5+2erQxE2S5O4OUlJM07xlSzPGLdvPP8OGDfDCC+DlZWF8wml5e5s8vn69GUFLvXrw3HNmTaItW6wOT9wEpS3aaisyMlLHy5Rnx/jXv0xv6YoVZtFuTCusa1dTkjl0CPz9LY5ROK30dDN6JigINm8Gj0sXzN6LLVqYnymp5zkVpdRmrXVkYedJy93VnT4Nb74JffvaEjvA8uWwapVplUliFwXx9jbVmK1bYdYsICTEzJVYuRIWLbI6PHGDpOXu6h55xIxl277dNL8wfatt25oVf/fuNftiC1GQzExT1UtLM/u7eOl003JPSzM/W9JCcBrSci8LVqwwU09ffNGW2MEs67t5s2mNSWIXReHpCa+/Dvv2ZW/O5O1txrz/+quZOyFcjrTcXVVysmlZKQXbtpmVoDBrtTdpAjVqwMaNZhNsIYpCa7jtNtNHs2ePqc4wciRMnmyWkGzZ0uoQBdJyd3//+AccPGjWA8lO7GD6Vk+cMCu4SmIXxaGUGeKemGjmMwFmOYLKlU2Sl+34XIr8+rui2Fgzi3D4cOjSxXZ4zx4zG/VPfzItMCGKq3VreOYZ+Phj82NGhQrw73+blrusO+NSpCzjai5eNHPGMzPN8IYKFQDzkbpLF1Oh2bMHqla1OE7hsi5dgmbNzM9QXBx4eWqzCPysWWbXpuhoq0Ms06Qs4460NtsoHT1qVgPLTuxg+r5WrzaNK0ns4mYEBZnG+tat8MEHmHrNJ5+YpS3uv980MITTk+TuSiZPhunTTb09V93l4EEzC7VHD7NrmhA36957oXdvswbd9u2Y3tVp00zD4oknstcqEM5MkruriI2Fp582E5XGjrUdTk83NXZvb/jyS5lMKBxDKfPzFBICgwebUVi0bw+vvWYaGLJypNOT5O4Kjh6Ffv2genXzi+XpaXvopZfMkMcvvjALQAnhKKGhZn7cnj3w7LPZB//6Vxg0CP7yF1i40NL4RMEkuTu7y5fN0gIpKeaXqUoV20Nz5sA778BTT8HAgRbGKNxWt26mNPPf/8LUqZjxtVOmQJs2pv6+bZvVIYrrkOTuzK5eNcXPHTvMkr633mp7aNs2swBkdLR8QhYla/x4s0XjyJFmlVH8/WH+fLPT9l13wYEDVoco8iHJ3VmlpprEvny5qbncdZftoePHTWdXcLBpvfv6WhincHteXjB7tpn13LevWa+IGjVg6VLT6dO1Kxw+bHWY4hqS3J1RWprZRWHxYvj0U7udlc6eNXn+4kWzYF/16hbGKcqMqlVNLvfwMD9/v/+O+SS5fDlcuWImWRw9anWYIhdJ7s7mwgWz2emCBfCf/8Bjj9ke+v136NzZfAqeN0+W+hClq0ED095ITIQ774RTpzA/hMuWwfnzZjRNQoLVYYpsktydydGjprj5889mJ5ynn7Y9dPy42dD4119Niz3X0u1ClJrISJg71zQwOnbMrsZERpqfWaXg9ttNa15YTpK7s1i5Etq1M1l8yRJ48EHbQ3v2QKdOcPKk+WjcrZuFcYoyr3t3+PFHUyLs0MH09xMebsbk1qsHvXrBhAlmYwFhGUnuVsvMhHHjTMYOCYF16+ya5bNnQ1SUKWuuWGFaS0JYrX1701gHk+BnzcJMtFizxszJePFFc3v2rKVxlmWS3K20Y4dpkr/6qhnXGBdnG+6Ynm7miQwcaA798otJ8kI4i1tvNUMjmzY185qeeQZSfYNNpv/gA1OeCQ83LRRZrqDUSXK3QlKS2dC6VSszruybb8zEkMBAwKw0EBn5xwSl1atl9qlwTrfcYlrwf/6z6f+PjobYOAWjRpkyTWioaaHcfbcMlyxlktxLU2qqadE0aGB21XjwQVNQf+ABwAw4eOYZ8wty9qzpuPrwQ9kqTzg3Hx8zkW7BAjOCJjraDPI6U7uV+TT6zjumptikCTz3XPYwG1HitNaW/GvTpo0uM86f13rCBK1vuUVr0DomRut162wPnzun9SuvaB0UpLVSWo8apfXFi9aFK8SNunhR6zFjtPb01DokxPxcnzmjtT52TOuRI80DAQFaP/OM1vv3Wx2uSwLidRFyrCT3kpKVpfXatVo//rj5YQatO3XSetky85jWetcurZ9/XuvgYPPwvfdqnZBgcdxCOMD27VoPGGB+rgMCtB49WuudO7XWe/dq/dBDWnt7m5ZMnz5az5ypdUqK1SG7jKImd9mJyZHS02H9ejMQfdYss9Own59ZM/XZZ6FVK06cMB9fp041A2O8veGee8ziTC1aWP0ChHCsnTvN1o/TppmBYW3bwrBhcHf0KWrO+9AsrXHypFlLo39/s75B9+5mxxCRr6LuxFSk5K6Uugv4N+AJfKG1/tc1j5cDvgLaAGeBwVrrwwVd0y2Se3Ky6f1cv95k6nXrzLoA3t5mOvYDD3Cp2wDWby/Pzz+biXw5L7lJE7OxxrBhsnOScH+nTpkGzeTJ2ePiMZNbe92VRaegrbRL+IyQpd+ZGdre3mZEQYcO5t9tt8kvSS4OS+5KKU9gH9AdOA7EAUO11rtynfMU0EJr/YRSaggwQGs9uKDrukxyz8qC06fN7NFff4Vdu8y/nTth/37bjvBXm7bkaHhvdtbrww7vluzY78uOHaa/NCvLLMEeFWWG/t59t9mjUjbWEGWN1uZXZ9Ei82/9evMrpBSEhWnCqyXSLDWBsNP/397dhch51XEc//5m9mWS3Uk222ytbnZNg14kjRaLREMFxfpSNdSbCFWUghe9sdBSS7AWRbzxoqAiChJUkFoQQcWglVix3hlJW1shbA2hqHlpbF528z5Jd/bvxZnNTDabdIw7c7LP/D5wmOc8c5j5c9j9z3nOc57n+TMbX/0d75idYpBLMD6e1l7ecQds2pQWJUxMpGVkPXbnvKVM7luBb0TExxv1xwEi4lstbXY32vxFUh9wFBiL63x4V5J7vQ6zs6mcP5+WIF6jzB6fofafU9ReP53K8bNcPHaaC0dPcfqNCtOs4SSjTGuU6ZENTK9ez2sDkxysv42DM1WOn7hy4dGGDbB5M9x5Z7oie+vWyysdzazhzJnmwe+ePWnctHDF5NjQOcYHjjFeP8j4uf2srR9lhJlmGYGR21Yw/NYqldGVVG4ZojJWpTJWZfDW1ZRGVqV/vmo1Lcpf5j8G7Sb3vjY+axw42FI/BLzvWm0iYlbSKeAW4Hh74bbvqa/t53tPXmQu1FJgLkrMkeoRStuUFpQB5qgwxxiBqFPmEgPUqFBvqyuAgP6zsKYf3lJNg4ctjQHE5GQakW/c6ERu1o5qNV2c3XpLjXPn0uUfU1Pp+cCHDw81ynr2HvkAJ05Avd5y2DvTKK8s/h0DXGSQi5SpUx7po9SfjqTL5XSXy8W2y+Urj6zntxfbdyPv79iRzrV1UjsZbbHJg4Uj8nbaIOlB4EGAycnJNr76aivWVLh19BylEleXMqgkSiVdrpdKSq/9ZUr9fai/j9JAswwOBZXhOpWRQSqrBqmsEJUKl8vgYDq3s2YNjI6m15UrPaVi1ilDQ3DXXalcTUSkH4CZmVSmp9Pr+fPp+Ta1GtTOzlKbvkBt5gK1M7PUzs5Sv/gGc+PD1EvpoH5uLr3Ol4X1efPzD63zEItt/y/vd+PalXaS+yFgoqW+DjhyjTaHGtMyq4GTCz8oInYCOyFNy9xIwNsfnWT7ozf2w2Bmy5+UjoyHh6935XYfUG2U3tTOFap7gXdKul3SAHA/sGtBm13AA43t7cCfrjffbmZmnfWmI/fGHPpDwG7SUsifRMQ+Sd8kLabfBfwYeErSAdKI/f5OBm1mZtfX1lnEiHgGeGbBvq+3bNeAzyxtaGZmdqN84zAzswJycjczKyAndzOzAnJyNzMrICd3M7MCcnI3MysgJ3czswJycjczKyAndzOzAnJyNzMroGzPUJV0DPhXB79iLR24n/wy5H5ocl80uS+alltfvD0ixt6sUbbk3mmSnm/naSVF535ocl80uS+aitoXnpYxMysgJ3czswIqcnLfmTuAm4T7ocl90eS+aCpkXxR2zt3MrJcVeeRuZtazCp/cJT0mKSStzR1LLpKelPSKpL9L+rWkkdwxdZukeyX9Q9IBSV/JHU8ukiYkPSdpStI+SQ/njiknSWVJf5P029yxLLVCJ3dJE8BHgX/njiWzZ4HNEfFuYD/weOZ4ukpSGfgB8AlgE/BZSZvyRpXNLPDliNgIvB/4Ug/3BcDDwFTuIDqh0Mkd+A6wA+jpEwsR8YeImG1U9wDrcsaTwRbgQES8GhGXgJ8Dn84cUxYR8VpEvNjYPkNKbON5o8pD0jrgU8CPcsfSCYVN7pLuAw5HxMu5Y7nJfBH4fe4gumwcONhSP0SPJrRWktYD7wH+mjeSbL5LGvzN5Q6kE/pyB/D/kPRH4LZF3noC+Crwse5GlM/1+iIiftNo8wTpsPzpbsZ2E9Ai+3r6aE7SMPBL4JGIOJ07nm6TtA14PSJekPSh3PF0wrJO7hHxkcX2S3oXcDvwsiRI0xAvStoSEUe7GGLXXKsv5kl6ANgG3BO9t/71EDDRUl8HHMkUS3aS+kmJ/emI+FXueDK5G7hP0ieBCrBK0s8i4vOZ41oyPbHOXdI/gfdGxHK6OdCSkXQv8G3ggxFxLHc83Sapj3Qi+R7gMLAX+FxE7MsaWAZKo52fAicj4pHc8dwMGiP3xyJiW+5YllJh59ztCt8HqsCzkl6S9MPcAXVT42TyQ8Bu0gnEX/RiYm+4G/gC8OHG38JLjdGrFUxPjNzNzHqNR+5mZgXk5G5mVkBO7mZmBeTkbmZWQE7uZmYF5ORuZlZATu5mZgXk5G5mVkD/BWdYzCGExe7pAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import time\n",
    "\n",
    "n_samples=200\n",
    "np.random.seed(42)\n",
    "start=time.time()\n",
    "U=np.random.rand(int(n_samples/2),2)\n",
    "R=np.sqrt(-2*np.log(U[:,0]))\n",
    "Box_Muller_sample=np.hstack((R*np.cos(2*np.pi*U[:,1]),R*np.sin(2*np.pi*U[:,1])))\n",
    "end=time.time()\n",
    "print(\"The time of sampling %d samples using Box-Muller algorithm is %.6f\"%(n_samples,end-start))\n",
    "\n",
    "start=time.time()\n",
    "gaussian_samples_randn=np.random.randn(n_samples)\n",
    "end=time.time()\n",
    "print(\"The time of sampling %d samples using np.random.randn is %.6f\"%(n_samples,end-start))\n",
    "plt.figure(figsize=(6,6))\n",
    "seaborn.distplot(Box_Muller_sample,\n",
    "                 hist=False,\n",
    "                 kde=False,\n",
    "                 fit=scipy.stats.norm,\n",
    "                 fit_kws={\"color\":'r'},\n",
    "                 label=\"Box-Muller algorithm\"\n",
    "                )\n",
    "seaborn.distplot(gaussian_samples_randn,\n",
    "                hist=False,\n",
    "                kde=False,\n",
    "                fit=scipy.stats.norm,\n",
    "                fit_kws={\"color\":'b'},\n",
    "                label=\"np.random.randn algorithm\")\n",
    "plt.legend(fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Rejection sampling\n",
    "In rejection sampling, we create a **proposal distribution** $q(x)$ which satisfies $Mq(x) \\geq \\tilde{p}(x)$ for some constant $M$, where $\\tilde{p}(x)$ is an unnormalized version of $p(x)$. We then sample $x \\sim q(x)$, which corresponds to picking a random x location, and then we sample $u \\sim U(0,1)$, which corresponds to picking a random height ($y$ location). \n",
    "If \n",
    "$$\n",
    "u >\\frac{\\tilde{p}(x)}{Mq(x)}\n",
    "$$\n",
    "we reject the sample,otherwise we accept it. The probability of acceptance is\n",
    "$$\n",
    "p(accept)=\\int \\frac{\\tilde{p}(x)}{Mq(x)}q(x)dx=\\frac{1}{M}\\int\\tilde{p}(x)dx\n",
    "$$\n",
    "So we need to select M as tight as possble.\n",
    "\n",
    "Suppose we want to draw samples from the posterior\n",
    "$$\n",
    "p(\\theta|D)=\\frac{p(D|\\theta)p(\\theta)}{p(D)}\n",
    "$$\n",
    "We can use rejection sampling with $\\tilde{p}(\\theta)=p(D|\\theta)p(\\theta)$ as the target distribution, $q(\\theta)=p(\\theta)$ as our proposal and $M=\\underset{\\theta}{max}\\,p(D|\\theta)$ is the MLE. We then accept points with probability\n",
    "$$\n",
    "\\frac{\\tilde{p}(\\theta)}{Mq(\\theta)}=\\frac{p(D|\\theta)}{p(D|\\hat{\\theta})}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "suppose we want to sample from a Beta distribution:\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$$Beta(x|\\alpha,\\beta)=\\frac{\\Gamma(\\alpha+\\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)}x^{\\alpha-1}(1-x)^{\\beta-1}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAFpCAYAAACf/JPiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmczdUfx/HXmc3MGPtgshNZSzKJSmXNEiLLkDBGYxQqa5YMiiyF7IQha4okW1l/lWwjJXt2StmNbcx2fn8cZJl97sx37r2f5+Mxj35z7/d+73v8+My553u+n6O01gghhHB8LlYHEEIIkTGk4AshhJOQgi+EEE5CCr4QQjgJKfhCCOEkpOALIYSTkIIvhBBOIsmCr5SapZQ6q5Tak8gxLymlflNK7VVK/c+2EYUQQtiCSurGK6XUC8A14AutdYV4ns8J/ALU01qfVErl01qfTZe0QgghUi3JEb7W+kfgYiKHtAGWaq1P3j5eir0QQmRCbjY4x2OAu1JqE5AN+Exr/UVSL/L19dXFihWzwdsLIYTz2Llz53mtdd7UvNYWBd8NqAzUAryALUqprVrrQw8eqJQKBoIBihQpQnh4uA3eXgghnIdS6kRqX2uLVTqngTVa6+ta6/PAj0DF+A7UWk/XWvtrrf3z5k3VLyghhBCpZIuC/y1QXSnlppTyBp4B9tvgvEIIIWwoySkdpdRC4CXAVyl1GggF3AG01lO11vuVUmuA3UAcMENrneASTiGEENZIsuBrrVsn45jRwGibJBJCCJEu5E5bIYRwElLwhRDCSUjBF0IIJyEFXwghnIQUfCGEcBJS8IUQwklIwRdCCCchBV8IIZyELZqnCSFui4mBY8dg/37zdfo0uLqCuzu4ucGjj8Izz0C5cuZxITKSFHwh0ujSJVi1Cr79FtasgatX/3suZ06IizO/CKKizH8BfHygWjXo1AmaNjW/EIRIb1LwhUgFrWH9epg4EVasgNhY8PODgAB49lkoWxbKlIEcOe5/zeHDsG2b+Vq5Elq1goIFISQEunSBPHms+5mE40tyi8P04u/vr6UfvrA3kZEQFgbjx8OBA+DrC4GB8Npr8PTT4JKCq2KxsbB6NUyYAD/8YM41bhy0aQNKpd/PIOybUmqn1to/Na+Vi7ZCJENkpBnNP/oovPWWmZKZMwdOnYJRo8y8fEqKPZg5/Fdege+/h99/h5IloW1bqF8fjh9Plx9DODkp+EIkIiYGpk0zhb5bN/PfDRtg+3Zo1w48PW3zPk88AT//bEb7mzdDhQqwfLltzi3EHVLwhYiH1uZCbMWKZn69WDEzZ/+//0GNGukz5eLqCl27wt69ZhXPq6/Cp5+aLELYghR8IR7w++9Qty40bGhW1ixZYkbfNWtmzNx6kSKwaZO5LtCrl/mFEx2d/u8rHJ+s0hHitr//hoEDYfZsyJXLXEDt0gU8PGz4Jpcuwa5d8O+/EBFh1nBGRUGBAlC4sPkqXhxvb3e+/NLk+fhjk23JEhtnEU5HCr5wetevw+jR5is6Gnr0gAEDTNFPs7NnzQL9tWshPNzclZWUnDmhQQNcmjRh+Pv1KFQoO2+/bS7oLlhgbuASIjXkr45wWrGxZqXNwIFw5gy0aAEjRkCJEmk88ZUrMHcufPWVmQuKizPzNFWqQHAwPPWU+T57dvPl5gZ//WWW/Bw/bi4UrFhhqruHB28FBxMZOoyeQ7Lj5WWWhaZ0RZAQAGitLfmqXLmyFsIKcXFar1ql9RNPaA1aV62q9ebNNjjxH39oHRKiddas5sTly2s9aJDWv/1m3jQlYmK0/vlnrYOCtHZ11drTUw99dpUG8xYpPZ1wHEC4TmXdlYIvnMqWLVq/8IL5m1+ihNaLFtmgeG7ZonW9euakWbJoHRiodXi4TfJqrbU+fFjrdu10nHLRfb0+06D1iBG2O72wL2kp+PLBUDiF/ftNz5pq1eDgQZg0yTzWqlUaVt5s22bukqpWDXbsgGHDTLe0WbOgcmXbhX/0UZgzB7X7dz4uOZNWLKLf+3GsXCZLd0TKSMEXDu3UKejY0dzItH49fPih6Wfz1ltpWPFy9Ci0bAlVq5pCP2KEmXvv39/0R0gvFSqgtm9jVsgOKrGL1q/dYv/6v9Lv/YTDkYIvHNKFC2YNe6lSMH8+vPOOqdMDB5q2CKly5Qr06WM6o61cCaGhZtVN375pOGkKeXriPeVTls28iLe+TuN6UVzcfjhj3lvYPSn4wqFcvw7Dh5uVNmPHQuvW8OefMGZMGgbfWsOiRVC6NHzyieludugQDB4M2bLZMn6yFe5Yh6WzrnAipiCvVz9B3K7fLckh7IsUfOEQYmJg+nQzoh8wAF56CXbvNksYixRJw4mPHIF69cxvjkKFTBOdsDDT09hiz3Z4jHGDLrEmqhZjn/0Ktm61OpLI5KTgC7umtbmv6fHHoXNnKF7cLH3/9lsoXz4NJ46KMh8VKlSALVtMP+Rt28A/VV1p002Xwflp+vJ1+kUOIrz2+6YRjxAJkIIv7NaWLVC9umkypjV8840p9s89l8YT//wzVKpkPio0bGiW83Trlin3JFQKZizIil8BFwIiw4io1xL++cfqWCKTkoIv7M6hQ9CsmdlZ6sgR0754zx5T+NPU3OziRbPnYPXqcO0afPcdfP11ppi+SUzu3DB/kRvHdDHeOvMBuuEr5mKGEA+Qgi/sxrVr0K+fmWVZu/a/JZbBwWnsL6O1aYVQpozpnNa7N+zbZ3YnsRPVq0NoqGJ+bACLdpU2F5ZjY62OJTIZKfgi09MaFi82qyFHjIDXXzeFfuBAyJo1jSc/dAhq1za7mZQoATt3mi2s0nzijNe/v9l5q6vXTP5Zvg0GDbI6kshkkiz4SqlZSqmzSqk9SRz3tFIqVinV3HbxhLM7fRoaNTJ3xPr6mt2gwsIgf/40nvjWLRgyxFzt3bkTpkyBX34xO57YKTc30wzuRpwnwUXWoId/bD4KCXFbckb4s4F6iR2glHIFRgLf2yCTEGgNM2aYlTYbNph19OHhZt4+zTZsMHsKDh5sLgYcOGB2GXGAFpSlS5sOD9+dfJK5BfqanspyEVfcluTfcK31j8DFJA7rBiwBztoilHBu//4LDRrAm2+axTK7d8N779lgkcyxY9C8OdSqZRbur1kDCxeCn59NcmcW77wDzz8P3a9+xF8R2cwcmMznC2wwh6+UKgg0BaamPY5wdj/8YGZVNm2CiRPNYLxkyTSe9Pp1M+FftiysXm2u9u7ZAy+/bIvImY6rq5n2io51pfNjG9AbNphts4TTs8Vn2HFAX611kkMIpVSwUipcKRV+7tw5G7y1cBTR0aYlzcsvQ548pifZ22+ncZZFa9NI5848R/PmplXmwIHg5WWz7JlRyZLmR165uwhfPfeZmb7audPqWMJiyrRXTuIgpYoBK7TWFeJ57hhwZ/WzL3ADCNZaL0vsnP7+/jo8PDyleYUDunDBNJ/csMEssRw7Fry903jSHTvM3MaWLebu2M8+s9EFAPsRG2tW7Zw+Fcd+VZ5c+T3Mn4tsjGvXlFI7tdapuuU7zSN8rXVxrXUxrXUx4GvgraSKvRB37N1rdv77+WezwmTatDQW+/37zV6FVaqY9pizZpmWCE5W7MFM7Xz+OZy/4ML7T642F0NkasepJWdZ5kJgC1BaKXVaKRWklApRSoWkfzzhyFasMC3lb9ww27i2a5eGkx09CoGB5q6sNWtM6+JDh8xjDrD6JrUqVYJ334Xp3xfj59qD4aOPTOEXTilZUzrpQaZ0nFtYmOliUKmSaXSW6u4F+/aZUevChWYheteu8P776bsRiZ25ft0sb/XOEsOuS8XJUjif+dSTptuThVUsndIRIqU+/dTsQlWrllmNk6piv2uXuQhboQIsXWqGsUePmn71UuzvkzUrTJ4M+w+5MarmGvj1V3Njg3A6UvBFhtHa9MLp1ctcpP3uu1RsFLV5s1mk/9RTsG6d6Wh54oQp9AUKpEtuR9Cggblb+aNvynOwZhcYOtTs/yicihR8kSG0NoV+xAjTt37BAsiSJQUv3rABatQwdxTd2TD8xAmzpl5G9Mkybpy5IB5y41N0TCz07Gl1JJHBpOCLDBEaamYRunUzbWuSddes1rBqlWlwX6uWWUM/Zsx/G4bnyJHesR2Kn5/pC7dpqxezGyyGr74yd7oJpyEXbUW6+/hjU5+Dgsw2hEkumomLM1dyP/rIzDcXKWIuxAYGgqdnhmR2VHFx8OKLsHev5kCOquTzuGxW7ST745awmly0FZnWhAmm2LdpY9bYJ1rstTYbjlSsaJqaRUSYdfSHD0OXLlLsbcDFxfzSvXZN0aP4N2bpqlzAdRpS8EW6WbbM3Oz66qvmpqpEp3F++cVM3bRoYW4RnT/f3EQVGAju7hmW2RmULWsuns/fWIB1z4Wa6yCnT1sdS2QAKfgiXfz6q2nS6O9vaneCS75PnIDXXjPF/vhx0xP5jz/MRwJZJ55u+vWDUqWgy18DiIx1N6udhMOTgi9s7q+/zKYlvr6wfHkCrRJiYsyC/HLlzJ2xQ4fCn3+aif5MuFm4o/H0NBfPDx93Z7j/UvjiC7PhgHBoUvCFTV2/bop9RIRZZx9vq/nwcDP079XLrL7Ztw8++MAutxW0Z7Vqmf1RRuyoyf5cz5plmhYt4hAZQwq+sBmtzRr7336DL780m0rdJzYWhg+HatXg3DlYssSsxila1JK8wnzI8vFRhPh+hf7xR3PhRTgsKfjCZmbMMPP1Q4aYOzvvc+qUGVIOGGDm7PfuNStxlIr3XCJj5Mtn1ub/+GcBZhcYAL17Q1SU1bFEOpGCL2zit9/MTVV16phlmPdZs8YstQwPh9mzTaOznDmtiCni0bGjuWbe+9ogzh+5DJMmWR1JpBMp+CLNIiLMaso8eWDevHuuuWpt5gwaNjQ3T+3aBe3by6g+k3FxMfdIXLnhQa8CC0zbiogIq2OJdCAFX6TJnXn7Y8dg0SIzRQBAZKRZQ9+rFzRtapqelSplaVaRsPLlzWzOnL/rsulCBfOLWjgcKfgiTRYvNoV+yBCoXv32gxcvmvn6OXPMXqqLF8sKHDswcCCUKAEh2eZz65MJcPas1ZGEjUnBF6n2zz/w1ltmN8G+fW8/eOaMadYSHm6W6oSGOvWOU/bE29v0zT94tSAjb3Y3K6qEQ5F/iSJV7kzl3LhhBvJubsCRI+bq37FjsHKlaXov7MrLL0NAAAx3GcChyevMndDCYUjBF6kyd665i3bYMChTBrPM8vnnzcW+DRugdm2rI4pUGjsWPLO60iVmAjp0sNVxhA1JwRcp9tdf0L27qe/vvINpclazpll989NPZo5H2C0/Pxgx0oUNugbzv4g1d0ILhyAFX6TYO++Ye3PCwsD1yKH/iv3GjaYVo7B7wcFQ1T+aHnzKxT4jrI4jbEQKvkiR1atNR4SBA6GkOmKKfWysmcYpXdrqeMJGXFxg2kx3Lqo89F1ZHbZtszqSsAEp+CLZbt6Erl1NXe/Z6rQp9pGRsH696XopHMoTT0CPbjHM4E1+6rJAGqs5ACn4ItlGjICjR2HyyKtkaVIPLl0ye6I+/rjV0UQ6CR3mQdHcEQTvCuHWqvVWxxFpJAVfJMuff5qC36ZVDDVH1zcPLFsGTz1ldTSRjrJmhSlhXhygLCPePCKjfDsnBV8kSWszlePpqfn0UpDZjnDePDOlIxxe/cbutK56lOFnOrB/3PdWxxFpIAVfJGn1ajNzM6T8V/j98AVMnGi6pQmnMW5pUbK6RNJ5gC9x0bFWxxGpJAVfJComxvQ/K+l7ibe2tDU9FN56y+pYIoPle8SVT0MO8dNNf2Z22WF1HJFKUvBFombMMPdVjTofhEeLV6W/ihPrML4yL2XdQe+wcpw5FWN1HJEKUvBFgiIiYFD/aF5w+YlXnzljmuZIIzSnpVxdmDbyMpFxHrzb/JTVcUQqyL9ekaAR/SM4d8mdT/OPRi3/Fry8rI4kLPbYW7X54JGZLN5enBXLZJRvb6Tgi3idPBTJ2MketHVbhP+6EffsbCKcmlL0nlSM8uzhraBIrl61OpBIiSQLvlJqllLqrFJqTwLPv66U2n376xelVEXbxxQZSmsG19+K1jBsmq/cRSvu4/FqAz4vO5bTF735oL+s2LEnyRnhzwbqJfL8MeBFrfUTwIfAdBvkEhY62Gcmc45W563nfqdIR2lzLB6gFNXGtKALUxg/yYXt260OJJIryYKvtf4RuJjI879orS/d/nYrUMhG2YQVVq0i9BMfvNyieX+JtDkWCXj5ZYY/vYxH1L8EvxlHdLTVgURy2HoOPwhYbeNzioxy8CC/txrOlwTw7nsu5MuvrE4kMiulyDG8LxPjuvD7bhfGjrU6kEgOmxV8pVQNTMHvm8gxwUqpcKVU+Llz52z11sIWLl+Gxo35IGogObPH0au/h9WJRGZXqxZNq1/gVc/VDB6sOXLE6kAiKTYp+EqpJ4AZQBOt9YWEjtNaT9da+2ut/fPmzWuLtxa2EBsLbdqw9UhevouqR+++LuTMaXUokekpBUOHMjGyE25xUXTpIr3VMrs0F3ylVBFgKfCG1vpQ2iOJDNe/P6xezcBSi8iXz2xfKESyvPQSBWuUZoRHKGvXmp56IvNKzrLMhcAWoLRS6rRSKkgpFaKUCrl9yCAgDzBZKfWbUio8HfMKW/vySxg1ip9f/YT1BwrRty/4+FgdStiVIUMIuTqKakX/pkcPuJjgEg9hNaUt+gzm7++vw8Pld4Ol9u6FZ56BihWp6/UTv//hwrFj4O1tdTBhd+rW5ffwaJ66soHgYMWUKVYHclxKqZ1aa//UvFbutHVWV65As2aQLRtb+i5j7XoXeveWYi9SacgQKl7aRLdqO5k2DWQslzlJwXdGcXHQoYPZr3DxYj6cmhdfXwgJSfKVQsSvWjWoX58he5uTP18cb71l1gKIzEUKvjMaOdJsT/jJJ+zwrM7q1dCzp8zdizQaMoQcl0/wSfXl7NhhWmuLzEXm8J3N2rVQrx60bAkLFtC4iWLzZjh+HLJlszqcsHuNG6N//IkaFc6xe58bBw+CrMC2LZnDF8lz4gS0bm2aoc2Ywa7fFN99Bz16SLEXNjJkCOrKZSZVnG72UxhkdSBxLyn4ziIyEl57DaKjYelSyJqVYcMgRw6zQbkQNlGpEjRtSvl5/Xg7KJLp0+GPP6wOJe6Qgu8sunaFnTth7lwoVYoDB0zd79rVFH0hbGbwYIiIIDTbGHLkgPfekztwMwsp+M5g9myYORMGDIDGjQEYNQo8PeGdd6yNJhzQE09A8+bknj6CIX2us349fPed1aEESMF3fHv3wltvQc2aMGQIACdPmoH+m2/KBTWRTgYNgqtXCYkYRZkyZhVYVJTVoYQUfEd2/bpZjZMtG8yfD66uAHz6qXm6Z08LswnH9vjj0KIF7hPHMmbIVQ4fhokTrQ4lpOA7sq5dYf9+U+z9/AA4dw4+/xzatoUiRSzOJxzboEFw7Rr1d4+kXj0YOlT67FhNCr6j+uILM3c/cCDU/m+bws8+Mwt2+ia4a4EQNlKhArRoAZ99xqj+l4mIMPf8CetIwXdE+/dDly7w4osQGnr34YgI87G6WTMoU8bCfMJ5DBoE16/z+JrRtG0L48fDX39ZHcp5ScF3NDdumFFV1qywYMHdeXuAqVNNz7R+/SzMJ5xL+fLmOtL48Qx55yKxsXfXDggLSMF3NN27w759ZieKAgXuPnzzJowZA3XrQuXKFuYTzuf2KL/4kk8ICYFZs+DgQatDOScp+I5k3jyz3r5fP1PZ7zF7Nvz7r4zuhQXKlYNWrWDCBAZ2uYCnJ3zwgdWhnJMUfEdx6JDpb1y9+kOfmWNizI1WVauaaX0hMtwHH8D16+T74hN69oSvvpKe+VaQgu8IoqPh9dchSxYzb+/mdt/TixaZbpj9+5t9p4XIcOXKQUAATJhAz/bnyZ3bdGAQGUsKviMYOtQMl6ZNg0KF7nsqLg5GjDAr5Bo2tCifEGBG+TdukH26GeWvXGnaO4mMIwXf3m3eDMOHmx2smjd/6OkVK0x3hfffBxf5f1tYqWxZ05574kS6BpwnVy4zVhEZR0qAPYuIgDfegKJFzR1VD9Aahg2D4sXNNTMhLPfBB3DzJtmnjea992D5cti1y+pQzkMKvj17912zqcncuZA9+0NPr18P27eb0f0D0/pCWKNMmbuj/G4B58iRAz780OpQzkMKvr1auRLCwkw1f+65eA8ZPtwsxW/fPoOzCZGYDz6AyEhyfj6ad9+Fb76B3butDuUcpODbo0uXIDjYXIlNYA+5LVtg40bo1css3hEi0yhdGtq0gUmTeKfNObJnl7n8jCIF3x716GHuogoLS7CaDxsGvr7m94IQmc7tUX6uz0fRrZvZfU3uvk1/UvDtzapV5rbZvn3BP/6N63/7zcz4vPuuaakjRKbz2GPm3pFJk+gecBYPj//2aRDpRwq+Pbl82QzZy5VLcCoH4OOPzTXct9/OwGxCpNTAgXDrFvlmjyIwEObMgX/+sTqUY5OCb0/efx/OnEl0KufgQXPb+ttvQ86cGZxPiJR47DGzE8/kyfR44xzR0TBhgtWhHJsUfHuxZYu5k7Z7d6hSJcHDRowwm5O/914GZhMitQYOhKgoSi0ZQbNmMHkyXL1qdSjHJQXfHkRHQ+fOpm1CIssZTpwwDTODg2VzcmEnSpUyo/wpU+jd8QKXL8OMGVaHclxS8O3BuHHwxx/m8262bAkeNmqUaY7Wq1cGZhMirW6P8p9ZN4wXXoCxY80YR9ieFPzM7sQJ01awcWN49dUED/vnH9MKv0OHh/qnCZG5lSxpWoRMmUKfNy9x6hR8+aXVoRxTkgVfKTVLKXVWKbUngeeVUmq8UuqwUmq3Uuop28d0Yt26mWF7Elezxowxo6I+fTIolxC2NHAgREdTP/xDypQxH2q1tjqU40nOCH82UC+R5+sDpW5/BQNT0h5LAGbN/XffmRF+kSIJHnbxIkyZYtqNlyyZcfGEsJlHH4V27XCZNoXu7a+wcyds3Wp1KMeTZMHXWv8IXEzkkCbAF9rYCuRUSj1iq4BOKyrKLLV57DGzMicRn30G167J9oXCzg0YANHRvHFqODlywPjxVgdyPLaYwy8InLrn+9O3H3uIUipYKRWulAo/d+6cDd7agU2YYLYtHDsWPDwSPOzSJfPxt1kz01pHCLv16KPQti0+YRMIan2Dr7+Gv/6yOpRjsUXBj2/TvHhn37TW07XW/lpr/7yybjBh//5rll82aGC+EjF2rGmLHxqaQdmESE/9+kFkJG8zidhYmDrV6kCOxRYF/zRQ+J7vCwF/2+C8zmvgQLhxw1yJTcTFi2Z037w5PPFEBmUTIj2VLg0tWlBi/oc0rh/FtGkQGWl1KMdhi4K/HGh3e7VOVeCK1vqMDc7rnH791ayvfOcd85c/EWPGmLsSZXQvHMqAAXD1Kt3zLuLcOVi0yOpAjkPpJNY+KaUWAi8BvsC/QCjgDqC1nqqUUsBEzEqeG0Cg1jo8qTf29/fX4eFJHuZctIZatcxNVocPQ44cCR564QIUKwb168PixRkXUYgM0aQJ+sefePyRc3h4urJzp1mdLEAptVNrHX+r3CQkufGd1rp1Es9rQPoy2sL335tdS8aPT7TYg2kle/26jO6FgxowALX8GbrVXE/I0rps2wZVq1odyv4lOcJPLzLCf0BcHFSqZNZX7t+f6Mqcc+egRAlo2FA+7goH9vLLXP31TwreOkLTpoo5c6wOlDmkZYQvrRUyiwULzMaew4YlWuzB9Lu/cUNG98LBDRhAtvPHaFvxD7780ixSEGkjBT8zuHXLrMx56ilo2TLRQ0+dMi1k27eHsmUzKJ8QVnjhBahenZBDPbl1y2z0JtJGCn5mMGWKaZI2ciS4JP5/yYcfmmu7MroXTmHgQJ44u45nH/2HqVOlv05aScG3WkQEfPQR1K0LtWsneuihQzBrFoSEQNGiGZRPCCvVqQNPP02XiFH8+adZ0yBSTwq+1SZMMGsshw1L8tDQULObVf/+GZBLiMxAKRg4kObnJpPHJ5Ip0poxTaTgW+nKFbO+slEj8E/8ovtvv5kVOe++C/nzZ1A+ITKDV17B8/HHCPSYz7JlmjNyW2eqScG30vjxpvvZ4MFJHtqvH+TKJbtZCSfk4gIDBtD54sfExChmzbI6kP2Sgm+Vy5fN6L5JE7M6JxHr1sGaNWYhT86cGZRPiMykeXNKPuZKzWzbmTlTExdndSD7JAXfKuPGmSmdJEb3cXHQu7dpo/C23M8snJWrK/TuTdDVcRw7pti0yepA9kkKvhUuXTJ9jZs2hSefTPTQ+fPN/P3w4ZAlSwblEyIzatuWpvm3kNPtKjNnWh3GPknBt8K4cWY5ZhKj+8hI0ziwcmVo1SpjogmRaXl64vVeCK/HzGHJ13FcumR1IPsjBT+jXb1qLta++mqSTezHjzd31o4eneT9WEI4h5AQgrJ+ya0oFxYssDqM/ZEyktGmTTMXbJPYgPb8eTON07Ah1KiRQdmEyOxy5KDS289SiV+ZOfmW1WnsjhT8jHTrltm1pGZNqFIl0UNDQ03jzJEjMyibEPbi3XcJcp3Drn1Z+PVXq8PYFyn4GWnOHDhzJsnR/Z49Zi/PkBAoXz6DsglhLx55hDZtIAuRzJx4w+o0dkX64WeUmBgoU8bcPbV9e4Lb92gNL78MO3aYTa/y5MngnELYgwMHeL3sTlZmeY0zlzzx8rI6UMaRfvj24Ouv4cgR0wgnkb3aVq2CtWvNAh4p9kIkoEwZgqru48otT5YuirI6jd2QEX5G0NrsZhUVZeZrElhyEx0NFSqY3wd//AHu7hmcUwg7Erd+IyVrF6VYaU82HChgdZwMIyP8zO6HH+D336FPn0TXV06aZFogf/qpFHshkuJS8yU6PrKajQcLcORP6bWQHFLwM8KYMfDII9CmTYKHnD1rpnHq1oUGDTIumhB2Syk6vP8ILsQS9sFRq9PYBSn46W3PHjPC79o10b1q+/WD69fhs88SneIXQtyjUMgr1MuyidnLchAba3WazE8KfnobOxa8vc0aywRs3252snr3XbOQRwiRTB4eBLWI4K9befl+8hGr02RO06EmAAAgAElEQVR6UvDT07//wrx50KED5M4d7yFxcdCtG/j5wQcfZGw8IRzBK5+8RF7OMfOTi1ZHyfSk4KenSZPM0pt33knwkDlzzAh/5EjInj0DswnhIDzy56Ldk7tZfvJJ/t1zzuo4mZoU/PRy8yZMnmy2L3zssXgPuXIF3n8fqlWDtm0zOJ8QDiRoWAlicGduz9+sjpKpScFPL3Pnms3Je/ZM8JAhQ+DcObOPuXTDFCL1yjYoTtUc+wnbUAR9S27ESoiUmfSgtVlu89RTUL16vIfs22cKfadOpt+9ECJtAtvFsi+mNDtGbbQ6SqYlBT89bNxoKnr37vGusdTaPOXjA8OGWZBPCAfUakg5vNRNZk+8ZnWUTEsKfnqYOBF8fRPcpuqbb2D9ehg6FPLmzeBsQjioHLlcaFb5JAvP1iTyx+1Wx8mUpODb2okT8O238Oab4On50NM3bkCPHvD449CliwX5hHBggYMKc5lcLOu3zeoomVKyCr5Sqp5S6qBS6rBS6v14ni+ilNqolNqllNqtlHLe5gBTppj/JnCj1ejR5nfChAng5paBuYRwAjUaelMk20XCtpQ2e0+I+yRZ8JVSrsAkoD5QDmitlCr3wGEDgcVa60pAADDZ1kHtws2bMGOG2a+2SJGHnj55EkaMgJYt4cUXLcgnhINzcYH27WCtrs2pTxdbHSfTSc4IvwpwWGt9VGsdBSwCmjxwjAbu3DaUA/jbdhHtyKJFZilmt27xPt2nj/nv6NEZmEkIJ9OhR240LnzxeaTZeEjclZyCXxA4dc/3p28/dq/BQFul1GlgFRB/xXNkWpt5mgoV4h2+//QTfPkl9O0b7+BfCGEjJUrAi+XPMzuiGXr5d1bHyVSSU/Dj69344K4prYHZWutCQANgrlLqoXMrpYKVUuFKqfBz5xzsFuht22DXLnj77YeWYsbGmmWYhQv/N8oXQqSfwB65OEwpNn/8o9VRMpXkFPzTQOF7vi/Ew1M2QcBiAK31FsAT8H3wRFrr6Vprf621f15HW484dSpkyxZvj4RZs+C338xUjre3BdmEcDLNW7ni43GLsPAK8OefVsfJNJJT8HcApZRSxZVSHpiLsssfOOYkUAtAKVUWU/AdbAifiIsXzXxN27bmbqp7XL4MAwbA88+bi7VCiPSXNSu0fC2WxbTk+oRZVsfJNJIs+FrrGKAr8D2wH7MaZ69SaqhSqvHtw3oCbyqlfgcWAh20VZvlWuGLLyAyEjp3fuipoUPh/HkYP142NhEiIwW+5c01svH1zCvmBhghm5inmdZQtqzpd//LL/c9deCAucEqMBCmT7conxBOSmt4rPANCv61nU2zjpl/iA5ANjG30v/+BwcPxju679HDzNl/9JEFuYRwckpBhxAv/sdLHB37rdVxMgUp+Gk1bRrkzPnQBP2qVbB6NYSGQr58FmUTwsm1a69QSjP7j6fMygknJwU/Lc6ehSVLzBaGXl53H46Kgvfeg9Klzd7lQghrFC4MdV6KYQ4diJs+w+o4lpOCnxZhYWYLw+Dg+x6eMAEOHTL7l3t4WJRNCAFAYLA7JynCxjknnf7irRT81IqLg88/hxdeMBdtb7twAT78EOrXN19CCGu9+irk9Ikm7EZLWOzc/XWk4KfW//4HR46YNsj3GDYMrl6FUaMsyiWEuI+nJ7Ru68YS1ZwrUxZYHcdSUvBTa8YMyJEDXnvt7kNHj5q9TwIDTUsdIUTmENhREak9+XJ7Mdi71+o4lpGCnxoXL5qLtW3b3nexdsAAcHc3N1sJITIPf38oVzqGMNXRTMU6KSn4qTF/Pty6ZXYgv23HDtMduWdPKFDAwmxCiIcoBYGd3Niqq7I/bKu5M94JScFPKa3NdE7lyvDkk3cf6t3brLfv3dvifEKIeLVtC64uccyOaApLl1odxxJS8FNq507Yvfu+0f3q1eYabmioaZgphMh8/PygQQPFXNcOxMz6wuo4lpCCn1IzZph5+9atAbM6c8AAs+nCAwt2hBCZTGBHxZnY/Pyw3hVOnUr6BQ5GCn5KXL8OCxaYNgo5cgDw9dfmju2hQ80FWyFE5tWwIfjmjiWMDjB3rtVxMpwU/JT4+muzyD4oCDDbZX7wgVmCGRBgcTYhRJI8PKBtO1eWqyZcmLnMXIBzIlLwU2LOHHj0UbObCWaAcOiQubPW1dXibEKIZAkMhCjtwYKjz8CWLVbHyVBS8JPr+HHYuBHatweluHULBg+Gp5+GJk2sDieESK4nnoCnnowlzCUIZs+2Ok6GkoKfXHfm+9q1A8yGJidPwvDhspOVEPamQ0dXdsU9ye8L9jpVQzUp+MmhtZnOqVEDihbl5k1T6F98EWrVsjqcECKl2rQBD/c4wq63gGXLrI6TYaTgJ8fmzaZRWvv2gFmZ+c8/ZkpHRvdC2J88eaBxY8V8lzeImjXP6jgZRgp+csyeDVmzwmuvERkJI0aYrsgvvWR1MCFEagV2VJyPy8PK9Z5mftYJSMFPyo0bpod28+bg40NYGPz9NwwaZHUwIURa1K0Lj+SLcao1+VLwk/LNN2btfYcO3LoFH38Mzz4LNWtaHUwIkRZubtAu0I1VNOSfGSucYk2+FPykzJkDRYvCCy8wZ465Gzs0VObuhXAEgYEQiyvzjj8Hv/xidZx0JwU/MadOwbp10L49UTEuDB8OzzwDdepYHUwIYQulS0O1KrGEqSB02Gyr46Q7KfiJmTfPfMxr144FC+DECTN3L6N7IRxHhyBX9umy7Fh42OHX5EvBT8idtffVqxNX/FFGjYKKFWVjciEcTatW4JUlltk3Wphrdg5MCn5Ctm2DgwehfXtWrID9+6FPHxndC+FocuSAZq+5sNDldSJnzrc6TrqSgp+QOXNM3/sWLRg5EooVM12RhRCOJ7Cj4nJcDpZtzOHQa/Kl4McnMtJsUNusGT/vzs4vv5i9at3crA4mhEgPNWpA0YLRDr8mXwp+fJYvh8uXoUMHRo4EX1/o2NHqUEKI9OLiAu2D3FlLHU7NXu+wa/Kl4Mdn9mwoVIg9eWuwYgV07w7e3laHEkKkpw4dQOPCF4ermW3sHJAU/AedOQPffw/t2jF6jCtZs8Lbb1sdSgiR3ooXh5eei2Y2HdDzHPPibbIKvlKqnlLqoFLqsFLq/QSOaamU2qeU2quUWmDbmBlo/nyIi+NMvUAWLjS7GebObXUoIURG6PCmO4cpxeY5hyE21uo4NpdkwVdKuQKTgPpAOaC1UqrcA8eUAvoBz2mtywPvpkPWjDF/PlSpwuQfShITY6ZzhBDOoXlz8PGMJuxCI9i0yeo4NpecEX4V4LDW+qjWOgpYBDy4qd+bwCSt9SUArfVZ28bMIHv2wG+/cbNle6ZONVsXPvqo1aGEEBkla1Zo2VKxmJZcn/2V1XFsLjkFvyBw6p7vT99+7F6PAY8ppTYrpbYqperZKmCGmj8fXF2ZR1vOn4d37fdzihAilQLfdOMa2fj6a+DmTavj2FRyCn5895Y+uGbJDSgFvAS0BmYopXI+dCKlgpVS4Uqp8HPnzqU0a/qKi4P589F1X2bcrOxUqmQ2ORFCOJfnnoNSBa8TFhkAK1ZYHcemklPwTwOF7/m+EPB3PMd8q7WO1lofAw5ifgHcR2s9XWvtr7X2z5s3b2ozp4+ffoJTp1hbsRf79sF770kbBSGckVLQIcSL//ESR6evszqOTSWn4O8ASimliiulPIAAYPkDxywDagAopXwxUzxHbRk03c2bB1mzMnZndfz8TEMlIYRzatfBBRcVx+wNReDiRavj2EySBV9rHQN0Bb4H9gOLtdZ7lVJDlVKNbx/2PXBBKbUP2Aj01lpfSK/QNhcZCV99xf5aXVmz1o233wYPD6tDCSGsUqgQ1Kl6jTlxbYn70nEu3iarO4zWehWw6oHHBt3zvzXQ4/aX/Vm5Eq5cYZJ+Cw8P6NzZ6kBCCKt16JaN1luys3Hyfmp1sTqNbcidtgDz5nE136N8sakwrVpBZru8IITIeK82VeT0vEnYHn+z+5EDkIJ/8SKsWsXcCiO4elVJGwUhBACentD6tWiW8BpXZi2xOo5NSMH/+mt0VBSTjjekcmWoUsXqQEKIzCLwnexE4sWXn0c4RAdNKfjz5vG/Iu3Yd9SLt9+WpZhCiP/4+0P5Ry4SduZl+P13q+OkmXMX/OPH4aefmJT9fXLnhoAAqwMJITITpSAwJAtbqcb+8WutjpNmzl3wFyzgLwrwzf4ydOxodjQUQoh7te2cFVcVy+yvvO2+g6bzFnytYe5cphf5iLg4RRcHWXYlhLCt/PmhYeV/mHutKTEbfrQ6Tpo4b8HftYuYA38yI6Il9epBiRJWBxJCZFaBvXw5QwF+GGXfO2E5b8GfN481rq/w9+WsBAdbHUYIkZk1aJoF3ywRhG0qbu7Mt1POWfBjY2HhQmbme998XGtodSAhRGbm4QFtG1zi25gGXFj4g9VxUs05C/6GDfzzj+a7f6vQrh24u1sdSAiR2XUcVIhoPJg/1j73dwJnLfjz5vGFZzCxcS4EBVkdRghhDx5/0hX//CeZuacK+uIlq+OkivMV/OvX0UuWMjPLWzz/PJQubXUgIYS9CArU7NZP8OuYTVZHSRXnK/jLl/Pz9Sc5dMWPTp2sDiOEsCcBfYrgqSKZGWaft+Q7X8GfN4+Z3t3Jlk3TvLnVYYQQ9iRnLkXzxw+x4O+XuPnnaavjpJhzFfyzZ7myZguLo5rQpo0ia1arAwkh7E3H3nm4Qk6WDrK/NfnOVfC//JJFcS24GeMhF2uFEKnyYpuClMhymlkr7G/jDOcq+PPmMcOrO48/brrgCSFESrm4QGDtU2y49gxHVx+0Ok6KOE/BP3SI3dtvEn6zPJ06SRtkIUTqdfioFIo4Zn9kX/P4zlPw589nJp3w8NC8/rrVYYQQ9qzQk7687LuTsG1liY2OszpOsjlHwdeayLlfMdetA82aKfLksTqQEMLeBbW+yenYAqwbv9fqKMnmHAV/yxaWHXuCSzHZ5WKtEMImGoU+RR7OM3NylNVRks05Cv7cucxwCaZY0Thq1rQ6jBDCEWTJ48Mbpbez7OjjnP/bPoq+4xf8W7c4tmAL6+Nq0jHIBRfH/4mFEBkkqFtW01At9JDVUZLF8cvfqlWERTRDKU2HDlaHEUI4kgrBz/K026/M/CobWludJmkOX/Bj58wjzCWIei9rChe2Oo0QwqG4uxP03EH+uFKUnT9etzpNkhy74F+4wA8rojgdV5CgTo79owohrBHQvwRe3GDW0My/Jt+xq+DixcyI7UDeXNE0amR1GCGEI8pRpwrNfdaw4MeC3LxpdZrEOXTBPzvzO5bTmHaBbnh4WJ1GCOGQlKLjq5e4EuPD0lmXrU6TKMct+H/+ydydZYnBnaBO0kdBCJF+Xuz/HI9ymJmfXbU6SqIctuDrufOYQSee9b9F2bJWpxFCODJVtgyBj3zPxj8Lc/So1WkS5pgFX2u2zNjLAcoS1CWL1WmEEE6gfSd3XIgl7NMLVkdJULIKvlKqnlLqoFLqsFLq/USOa66U0kopa5sPb97MjDMN8PGMpmVLS5MIIZxEoZBXeJnvmT3PjdhYq9PEL8mCr5RyBSYB9YFyQGulVLl4jssGdAe22TpkSkXM/IovaUVAK/DxsTqNEMIpFChA0OM7OB2Rg7U/ZM67sJIzwq8CHNZaH9VaRwGLgCbxHPchMAqItGG+lIuM5Msv4QZZCQpxtzSKEMK5NOpWDF/OMfOTS1ZHiVdyCn5B4NQ935++/dhdSqlKQGGt9YrETqSUClZKhSulws+dO5fisMmyciUzb7amfNFrPPNM+ryFEELEx6Plq7zhuoBvN2Xn/Hmr0zwsOQU/vjWNdz+vKKVcgLFAz6ROpLWerrX211r7582bPvtB7pm4iW1UJaibt+xqJYTIWDly0LHmCaLj3Jg3J/NN5Cen4J8G7u1CUwj4+57vswEVgE1KqeNAVWC5JRduz59n5v9K4u4SwxvtHXMBkhAic6vw9otUYRuzJl7PdA3VklMVdwCllFLFlVIeQACw/M6TWusrWmtfrXUxrXUxYCvQWGsdni6JE3Fr3lfM1a/zau1r+Ppm9LsLIQRQvz4dsy7mj+PZCc/wKpi4JAu+1joG6Ap8D+wHFmut9yqlhiqlGqd3wJT4duIpLuBLUI+cVkcRQjgrDw8C2rqZhmpTblmd5j5KW/SZw9/fX4fb8tffwYO8XOY4B3JW4+j57Li62u7UQgiRIr/+SrvKe/jWsxVnLmTB29t2p1ZK7dRap2rK3GEmuk+M+4a11CEwyEWKvRDCWpUqEVR8IxGRWVi61Oow/3GMgh8TQ9hcNwACu8mdVkIIiynFC29VMA3VJmSejVEcouDHrvqesOstqPPkOYoWtTqNEEKAavs6HdVsNm3PypEjVqcxHKLgrxv1KycpSlDvPFZHEUIIw8+P9jVPmYZqM+OsTgM4QsE/f56Zv5Qlj+c1mrzmZnUaIYS4q2BII+qxhtmfR2WKhmp2X/DPT1vCMt2YN5pHkkU6IQshMpNGjeiYdTF/nffk+++tDmPvBV9r5k6KIBoPgvrKnVZCiEwmSxYatctFPv7l80lRVqex74Kvf93FzDP1eabYv1SoYHUaIYR4mEfnQAIJ47s1bvz9d9LHpye7LvjbPt7AXioQ9I4sxRRCZFIVK9KpwjZi41wIm2Vtcx37Lfg3bjBzeV6yukUSEJTV6jRCCJGgkl3rUYt1fD45ijgLF+zYbcG/NvcbFkU3o2Wdy2TLZnUaIYRIROvWBHvM4cSZLKxda10Muy34X44+yTWyEdQ/v9VRhBAicdmz82qAJ3nVOaZPjrEshn0W/D17mH6kJuXyn+fZ52SXEyFE5ufROZAOOozlK104c8aaDHZZ8H8bvortPEPn7p6yq5UQwj5Uq0anEhuJiXVh9mxrIthfwb95k2lLfPF0jeKNLrI6RwhhJ5TisbfrUIMNTJ9kzZ23dlfwr81bxvyo5rSsdZFcuaxOI4QQKdCuHV3cZnD8Lw9L7ry1u4K/cNQprpKdzoPkYq0Qws74+vJqgCd+6h8mj4/O8Le3r4K/dy/TD9eggt85qj0rk/dCCPvj3i2EYD2NVT+4cexYxr63XRX8XwcvJ5yn6fyOl1ysFULYpypVePOJ7bjoWKZNzdg7b+2nn/CVK0xblh8v11u0DbG/i7URERGcPXuW6OiM/xgnhKNyd3cnX758ZM+e3eooKVLo3eY07ricmdNeYfAQDzw9M+Z97abgX5k8n/kx7WjV6Bo5c9pXH+SIiAj+/fdfChYsiJeXF0o+ngiRZlprbt68yV9//QVgX0U/IIC33nmdb6404+uvoW3bjHlb+5jSiYtj9piLXMeHrqH2t6vV2bNnKViwIN7e3lLshbARpRTe3t4ULFiQs2fPWh0nZby8qBlcksc4yJTPMq5tsl0U/LjV3zPpfEuqlTpH5cpWp0m56OhovLy8rI4hhEPy8vKyy6lSly6d6cJUfgn34LffMug9M+Zt0uaHwb/wJ4/RdaD9LryXkb0Q6cNu/209+igd6vxFVnWd8WMz5i6szF/wDx9mYvgz5M96jeYBdnPJQQghkpSzTzDt9WzmL4CMmJXK9AX/yLBFrKIBnYM1Hh5Wp3FegwcPRil198vPz49XXnmF3bt3p+p8SikmTpxo04yHDh1i8ODBXL582abntcqmTZtQSrFnzx6ro4j0UqsW3R/7nqgYV6ZOSf8lmpm74F++zOT5OXBVcXTuJU3vrZYjRw62bNnCli1bGDduHIcOHaJOnTpcvHgxxefasmULLVq0sGm+Q4cOMWTIEIcp+MIJKEXpAc2pzyqmfBbFrVvp+3aZuuBfnzCLWdFtaV73KgUKWJ1GuLm5UbVqVapWrUpAQABffPEFZ8+eZc2aNSk+V9WqVcmfX9pjCEFAAO/mnss/l7KweHH6vlXmLfi3bjH3k3+5TC66fmC/F2sdWcWKFQE4derUfY9fvHiRzp07kz9/fjw9PXn22WfZtm3bfcfEN6Xz7bff4u/vj6enJ35+fvTp0+eh1Re7d++mUaNG5MyZEx8fH6pUqcLatWvZtGkTjRo1AqB48eIopShWrFiC2ffu3Uu9evXInTs3WbNmpWzZskyaNOnu8ytXrqROnTp3b+qpWrUqP/zww33nGDx4ML6+vmzbtg1/f3+8vLx4/vnnOXbsGGfPnuXVV1/Fx8eHsmXLsmHDhvteW6xYMXr16sWHH36In58fPj4+vP7661y5ciWRP3GIi4tjxIgRlCxZkixZsvDYY48xZ86cRF8jMjkPD+r0qkhZ9jHu4xvodJzZybQFP3buAj6N6MTTpa/w7LNWpxHxOXnyJGAK7B23bt2idu3arF27ltGjR7Ns2TLy5s1L7dq1+eeffxI81+LFi2nWrBlVqlRh+fLlhIaGMn36dPr163f3mAMHDvDcc89x5swZpk6dyjfffEPTpk05deoUTz31FJ988gkAS5cuZcuWLXzzzTcJvl/jxo1xdXVl3rx5LF++nG7dunH16tW7zx87doxGjRoxd+5clixZwrPPPkv9+vXZvHnzfee5ceMGwcHBvPfeeyxcuJCTJ0/yxhtv0Lp1a55//nmWLl1KwYIFadGiBTdu3LjvtQsXLmTdunV8/vnnjBkzhpUrV9KpU6dE/8y7devGRx99RHBwMCtXrqRp06Z07NiRFStWJPo6kbmpkM684zGVX/d788BfMdvSWlvyVblyZZ2g2Fi9pGA3DVov/jIu4ePsxL59+6yOkGahoaE6T548Ojo6WkdHR+vDhw/r2rVr6yeffFJHRkbePW7GjBna3d1dHzp06O5j0dHRukSJErpXr153HwP0hAkTtNZax8XF6SJFiugOHTrc954zZ87Unp6e+vz581prrQMCAnTBggX1jRs34s343XffaUAfO3Ys0Z/l3LlzGtC7d+9O1s8eGxuro6Ojdd26dXVgYOB9fyaA3rRp093HJk2apAE9ZMiQu4/t3btXA3rVqlV3HytatKjOlSuXvnr16t3H5s2bp5VSd/++bNy4UQP6jz/+0Fpr/eeff2qllJ49e/Z9+d544w3t7++frJ/FUTnCv7HrIT10Li7oZg3i//t9BxCuU1l3k7XOUSlVD/gMcAVmaK1HPPB8D6ATEAOcAzpqrU+k+pfQylWM/qs1JfJdpdlrDnqx9t13ybC7LR705JMwblyKX3bhwgXc3d3vfp8nTx527NhBliz/tbpYt24dlStXpnjx4sTE/Ld354svvkh4eHi85z106BAnT56kZcuW972mZs2aREZGsmfPHl588UU2bNhA27Zt03wTW+7cuSlcuDAhISF0796dGjVqkC9fvvuOOX36NAMGDGDdunWcOXMGfftz9nPPPXffcR4eHlSvXv3u9yVLlryb/cHH7rQAuKNOnTr4+PzXF6pZs2a0bduWHTt2ULZs2Ydyr1+/HhcXF5o2bXrfn1OtWrVYuHAhsbGxuLq6pujPQmQe3r3fJmTadEas6sOff0KpUrZ/jySndJRSrsAkoD5QDmitlCr3wGG7AH+t9RPA18CotITaPHA1W6lGjwHeyN/fzCNHjhzs2LGDrVu3Mm3aNKKiomjTpg1xcXF3jzl//jxbt27F3d39vq+wsLCH5vrvfQ1AgwYN7nvNnamiO6+7cOECjzzySJp/DhcXF3744Qf8/Pzo2LEjfn5+VK9enV27dgFmnrxx48b88ssvDB06lI0bN7Jjxw7q169PZGTkfefKli0bLi7//TPyuL12OGfOnA899uBrH/wl4+XlhY+PD2cS2PD0/PnzxMbGkiNHjvv+nDp06EBMTEyCrxN2okQJ3mn+Fx5EMXrozXR5i+SM8KsAh7XWRwGUUouAJsC+OwdorTfec/xWIPWtgLZuZfTuuuTJepPATg7cjiAVI2yrubm54e/vD8AzzzyDl5cX7dq146uvvqJVq1aAGT37+/szZcqUh15/7yeBe+XOnRuA6dOnU6lSpYeev1P48+TJY7OiVqZMGZYsWUJ0dDQ//fQTffv2pWHDhpw+fZrDhw+za9cuVq9eTb169e6+5uZN2/4jfLD/y82bN7l27VqCv9Ry586Nm5sbmzdvvu+XzB0P/gIR9if/h13p+FUYMxe+yeCR2Hx1YnIu2hYE7h2anb79WEKCgNWpDXSg1wyW04S3u7ni7Z3as4iM0LZtW8qXL8/IkSPvPlarVi0OHz5MkSJF8Pf3v+/r8ccfj/c8pUuXpmDBghw/fvyh1/j7+5MnT5675168ePFDI+U7EhpJJ8bd3Z2aNWvSo0cPzpw5w+XLl+8W9nt/QZ04ceKhC7ZptXbtWq5du3b3+6VLl6KUuvtL9UE1a9YkNjaWK1euxPvn5CF3Jtq/0qXp1fggMbGKscNvJH18CiVnhB9fo4p4Fw4ppdoC/sCLCTwfDAQDFClS5OEDtmzhk81V8XSLpmsP+cub2Sml6N+/P6+//jrr16+nVq1atGvXjqlTp/LSSy/Rq1cvSpQowYULF9i+fTt+fn689957D53HxcWFTz/9lDfeeIOIiAjq16+Ph4cHR48eZdmyZXz99dd4e3sTGhrK008/zQsvvEDPnj3JkycPu3btIk+ePHTs2JHSpUsDMG3aNAICAvD29o73l8zu3bvp1asXrVq1okSJEly6dImRI0dSsWLFu8s0CxUqRM+ePfnwww+5evUqoaGhFCyY2Dgn5by8vGjYsCG9e/fmzJkz9O7dm6ZNm1Ku3IMzpkbp0qUJCQkhICCAPn364O/vT2RkJHv37uXQoUPMmDHDpvmENUp8HEzA8kVM/bw5/T/Etnt3J3VVF6gGfH/P9/2AfvEcVxvYD+RLztXi+FbpnKzeRnsQqUOColJ0dTuzc4QVBHdW6TwoJiZGlypVStetW/fuY5cvX9bdu3fXhQoV0u7u7rpgwYK6adOm+ueff757DKAnTpx437lWrVqln3/+ee3t7ashsn0AABPESURBVK2zZcumK1asqAcMGKCjo6PvHvP777/r+vXrax8fH+3j46OrVKmi161bd/f5Tz75RBcpUkS7urrqokWLxvuz/Pvvv7pt27a6ePHiOkuWLDp//vw6ICBAnzhx4u4x27dv108//bT29PTUJUuW1GFhYbp9+/b63r+38f2ZPLiy5t6f986qJK3NKp0ePXro0NBQnS9fPu3t7a0DAgL0pUuXEj1XXFycHjt2rC5Xrpz28PDQvr6++oUXXtBz5syJ92d1Fo7wb+xev9ftpUHrj/pff+g50rBKJzkF3w04ChQHPIDfgfIPHFMJOAKUSu4bP1Twf/5ZhzBZu7vG6OPHbfFHlnk42l/GtIqIiNCAXrBggdVRLFO0aFHds2dPq2M4DIf7N7Z7t27ACp3X+6q+/kDNT0vBT3IOX2sdA3QFvr89gl+std6rlBqqlGp8+7DRgA/wlVLqN6XU8pR+0jjZdxIzCaJj+ziKFk3pq4W9OHDgAB9//DFKKZ566imr4wiROT3+OO+/sIVzN3yY9um1pI9PpmStw9darwJWPfDYoHv+d+00pfjpJ4ZvfgFcXekfKuswHVm/fv3Yvn07o0ePvjvnLoR4WPVpbalVdj0jPn6aN98DHxts5W19g3mtOdFnErP4gk4d4yhSRAq+I0us3YEzOX78uNURRGZXpgwfNvmWZ7+txcQh53l/tG+aT2l9L52VKxm2tSbK1ZX+g9yTPl4IIZxEtantaei6mlGfeZJEX71ksbbgR0dz/J2xhBHIm8FQqJClaYQQInPx82No0EkuRfswtufpNJ/O2oI/dSofHO2Aq7sL7/eXqRwhhHjQU5++zmtZVjBmdi4unE9b72TrCn5sLNsHLmceb9Cjp4uM7oUQIj4+Pgzpe4NrsV6MevNQmk5lXcH/+ww9IkLJlzuafv3tdNd5IYTIAOU/aEabnKsY/23a1qxbVvAvnY1mM8/z0Qh3sjloB2QhhLAJNzc+nuWH0nFJH5sIywr+aQryRLloOna0KoFIqdmzZ1O5cmWyZctGrly5qFSpEj169EjROe5s7Sf+E992j5nRqFGj2LRp00OPp3f+5s2b89JLL6Xb+e1F4ab+vF95bZrOYVnBj8KDMePdpd+9nfj444/p1KkTL7/8MkuXLuWLL76gSZMmLF+espuqv/nmG7p3755OKUV6Sqjgi4zT+9vqSR+UCMtuvMqRA2rVsurdRUpNnDiRzp07M3z48LuPNWrUiNDQ0BSdJ75+90KI5PEqmDtNr7dshF+4sFXvLFLj8uXL+Pn5PfS4UvdfcL958yZ9+vShaNGiZMmSheLFi9+3EfmDUzodOnTA39+fZcuWUaZMGTw9PXn++efZt+/u/jq0aNGCGjVqPPTeoaGh5M+fn+jo6AQzd+rUiQIFCuDp6UmRIkV488037z5/4MABAgICKFy4MN7e3pQvX55x48bdt4PXpk2bUEqxfv16mjRpQtasWSlVqhQ//PADsbGx9O7dG19fXwoWLMiYMWPue//k/GwJ+fbbb/H398fT0xM/Pz/69OmT4M95x//bu/egqup2gePfBwLEEAxlYxIKZlGIqEVYZgcL0cTyMmmlvh5Q8lIenffYHHA8Sb456XhpdEz01bIotfSUeSlNisQ8lYod0saytzqGmPfLm9oRLeE5f2xYsJHLQmGj8PvM7Jm91v7ttZ71sNezF2ut/fv16tWLIUOG8OabbxIeHo6fnx8jR47k0qVL5ObmEhsbi5+fH7169bIGoC918eJFUlNTCQ0NxcfHhy5durB5c1lvKmFhYZw+fZq//e1viAgi4nK0X1RUxNSpUwkKCsLhcDBhwgQuXbrkso49e/YQHx9P8+bNueWWWxgxYgTHjx93aXPo0CESExPx9fUlLCzMdPlc166217VrfVQ7iHkj0xh68uvZs6cGBQVpZmamNah4RcXFxZqQkKB+fn46e/Zszc7O1rfeekufeeYZq03FXiKTkpK0devWGh4eritXrtS1a9dqVFSU3nbbbVpYWKiqqlu2bFER0QMHDrisKywsTCdPnlxlzKNGjdKIiAhdvXq1btu2TVesWKFjxoyxXs/Oztb09HTduHGj5uTk6Pz589Xf319nzpxptSntovj222/XOXPmaFZWlsbHx2uLFi10/PjxOm7cOM3KytJJkyYpoDt27KjVtqle2XXymjVr1MPDQ5999lnNysrSxYsXa0BAQI29a8bFxWlISIjGxcXphx9+qBkZGert7a1jxozR6OhoXblypa5bt05DQ0O1b9++Lu/t37+/BgUF6eLFizUrK0tTUlLU09NTv/nmG1VVzcvL04CAAE1JSdEdO3bojh079OzZs1b8oaGhmpSUpFu2bNE5c+aop6enzp4921r+iRMnNCAgQO+//35dt26drlixQkNCQrRz58566dIl62/arVs3DQ0N1VWrVln5atu2rcbFxVW77Y1hH7OL+uweub4epuDfWPbu3avh4eEKqIhoZGSkTps2zdrpVZ2FGdANGzZUuZzKCj6gX375pTUvPz9fPT09dcmSJaqqWlRUpO3atdP09HSrzWeffVZpv/PlderUSRcuXGhr+4qLi/XPP//Ul19+WcPDw635pQV/+vTp1rzvvvtOAX344YeteUVFRRocHKypqam12jZV14JfXFys7dq10+TkZJf4li9frs2aNavyy1bVWfADAgL0t99+s+YNHTpUAf3888+teRkZGQro/5X0u5udna2Abtu2zWV5Dz30kA4ZMsSabtWqlb744otXrBfQhx56yGXewIEDtXv37tZ0WlqaBgQEuHxedu3a5dJN9qZNmxTQnTt3Wm1K82UKfplrKfgN33laE/XXv8KePQ2z7q5daz+kbnR0NPv37+eTTz4hKyuLrVu3MmPGDFavXk1eXh5+fn5s3bqVwMBABgwYUPMCy3E4HPTo0cOabt++Pffeey+5ubmMHz8eDw8PkpOTefvtt5k+fToiQmZmJjExMURFRVWznV2ZO3cunp6e9O7dmzvvvNPl9YsXLzJr1ixWrVpFQUGByymTy5cvc9NNZbtHfLkLTh07dgScQw6W8vDwoEOHDhw+fLhW21bRjz/+SEFBAU8++SSXL1+25j/yyCNcvHiRffv2ERdX6YByAMTExBAQEOASq7e3Nz179rwi/iNHjtCxY0eys7Np06YNDz74oMs64+PjyczMrHJd5fXp08dlOjIykq+//tqazs3NpU+fPvj7+1vzYmNjCQsL44svvmDYsGHk5uYSHBxM9+7drTal+TLqRsN3nmbcMHx8fHj88cdZtGgR33//Pa+//jo//fQTy5cvB+D06dNVDsBdncoG33Y4HC4Dlo8aNYqDBw+Sk5PD+fPnWbt2LaNruKd30aJFDBo0iJdeeomIiAjuuOMOVq9ebb2elpbGvHnzGDt2LJs3b2b37t288MILwJXj4rZs2dJ6Xjp2bPl5pfMrvs/OtpV36tQpABITE/Hy8rIepQO5Hzp0qNL3VRZnaUwtWrRwGfS84ti/p06d4tixYy7r8/LyYvr06TWur7r1ls/F0aNHCQ4OvuJ9wcHBnDlzBoBjx45VmS+jbpgj/AZS2yPs61FKSgqpqan88MMPALRq1arKQladEydOVDqvU6dO1nRYWBi9e/cmMzOTX375heLiYoYNG1btclu2bMnChQtZuHAh3377LXPmzGHEiBFER0cTGRnJe++9x8SJE0lNTbXes2nTplrHXx0721ZeYKDzLoxly5ZVekdTaeGvS4GBgYSEhLB+/fo6X3apW2+9tdJcHD9+3DqCb9OmTZX58vX1rbfYmhJzhG/YUtmOePLkSc6ePWsducXHx3PmzBk++uijWi/7q6++sqYLCgrIy8sjNjbWpV1KSgpr165l8eLFDBo06IqjyupER0czd+5ciouLrS+owsJCfHx8rDZFRUUu/wHUBbvbVioiIoKQkBDy8/OJiYm54tGqVas6jQ+cf7djx47h5+dX6TpLVfYfjF3du3cnKyuL8+fPW/N2795Nfn6+dbrpvvvu4/jx4+zatctqU5ovo26YI3zDls6dOzNw4ED69OmDw+Hg4MGDzJs3j+bNm5OUlARAQkICffv2Zfjw4aSnp3PPPfdw9OhRtm/fztKlS6tcduvWrRk5ciQzZszA19eX9PR0HA4HycnJLu0GDRrEc889R15eHrNmzaox5p49ezJ48GCioqIQEV577TVuvvlmq9gmJCSQkZFBx44dCQwMJCMj44pbCa+V3W0r5eHhwSuvvMLIkSM5d+4c/fr1w9vbmwMHDrB+/Xref/99mjdvXqcxlv7dEhISSEtLo1OnTpw7d449e/ZY1zkA7rrrLjZt2sSjjz6Kn58fERERtLDZL8rkyZNZsmQJffv2JS0tjd9//50pU6bQuXNnnnjiCcB5GqtLly4MHTqU2bNn06xZMytfRt0wBd+wJT09nQ0bNjBp0iTOnDlDmzZt6NGjB2vWrLFOM4gI69atY9q0aSxYsICTJ0/Stm1bhg8fXu2y27dvz9SpU5kyZQoHDx4kJiaGd999l2bNmrm08/HxoV+/fmzfvp3evWseVfOBBx4gMzOT/Px8PD096datGx9//DG3lXTN+uqrrzJ+/HgmTJiAr68vSUlJDB48mLFjx15llq5+28p76qmn8Pf3Z+bMmbzxxht4enrSoUMHHnvsMev8e10SET744ANmzpzJggULKCgoIDAwkK5duzJx4kSr3dy5c5kwYQL9+/fnwoUL5OTk2O7yICgoiJycHJ5//nmGDRuGt7c3iYmJzJ8/39omEWHjxo2MHTuW0aNH43A4mDp1Kp9++ql1bcO4NuK8y8f9YmJitPxV/MZs//793H333Q0dxnUpOTmZffv2YeezcPnyZdq3b8/o0aOZMWOGG6K7NrXZNuPaNKV9TET+R1Vjam55JXOEb1z3/vjjD/bu3cs777zD6dOnGTduXEOHZBg3JFPwjevekSNHiI2NxeFwsHTpUuuUjGEYtWMKvtGg7PywJywsjIY69Xgt7P5oyTDcxdyWaRiG0USYgm8YhtFEmILvJjfiKQnDuBGYfcs+U/DdwMvLi8LCwoYOwzAapcLCQry8vBo6jBuCKfhu4HA4OHz4MBcuXDBHI4ZRR1SVCxcucPjwYfNrXJvMXTpuUNol7JEjR2octcgwDPu8vLwIDg526XbZqJop+G7i7+9vPpSGYTQoc0rHMAyjiTAF3zAMo4mwVfBF5FER+YeI/CwiUyp53UdE1pS8vktEwuo6UMMwDOPa1FjwRcQTyAD6AZHAMBGJrNAsBfinqnYE5gOz6zpQwzAM49rYOcKPBX5W1QOq+gewGhhYoc1A4K2S5+8D8SIidRemYRiGca3sFPwQoPxIxr+WzKu0japeBs4CdT8Wm2EYhnHV7NyWWdmResVfD9lpg4iMBUqHE7okIvtsrL8paA2YIX2cTC7KmFyUMbkoE3G1b7RT8H8FQstN3wYcqaLNryJyExAAnKm4IFVdBiwDEJGvr3bUlsbG5KKMyUUZk4syJhdlROSqh1Czc0pnN3CHiISLiDfwNLCxQpuNQFLJ8yHAVjV9CBiGYVxXajzCV9XLIvJvQBbgCbyhqt+JyEvA16q6EVgOrBCRn3Ee2T9dn0EbhmEYtWerawVV3QxsrjAvvdzzi8DQWq57WS3bN2YmF2VMLsqYXJQxuShz1bkQc+bFMAyjaTBdKxiGYTQR9V7wTbcMZWzkYrKIfC8i34rIZyLSviHidIeaclGu3RARURFptHdo2MmFiDxZ8tn4TkTecXeM7mJjH2knIjki8k3JfpLYEHHWNxF5Q0ROVHXrujgtLMnTtyJyj60Fq2q9PXBe5P1foAPgDewFIiu0eQ74e8nzp4E19RlTQz1s5uJhoHnJ82ebci5K2rUAtgM7gZiGjrsBPxd3AN8At5RMOxo67gbMxTLg2ZLnkUB+Q8ddT7n4F+AeYF8VrycCH+P8DdT9wC47y63vI3zTLUOZGnOhqjmqeqFkcifO3zw0RnY+FwAzgDnARXcG52Z2cjEGyFDVfwKo6gk3x+gudnKhQOnAEgFc+ZugRkFVt1PJb5nKGQi8rU47gZYicmtNy63vgm+6ZShjJxflpeD8Bm+MasyFiHQDQlX1I3cG1gDsfC7uBO4UkS9FZKeIPOq26NzLTi6mA38RkV9x3jk40T2hXXdqW0+A+h/xqs66ZWgEbG+niPwFiAHi6jWihlNtLkTEA2evq8nuCqgB2flc3ITztE4vnP/1/beIRKnqb/Ucm7vZycUwIFNVXxGRB3D+/idKVYvrP7zrylXVzfo+wq9NtwxU1y1DI2AnF4hIb+A/gQGqeslNsblbTbloAUQB20QkH+c5yo2N9MKt3X1kg6r+qaq/AP/A+QXQ2NjJRQrwXwCqugNohrOfnabGVj2pqL4LvumWoUyNuSg5jbEUZ7FvrOdpoYZcqOpZVW2tqmGqGobzesYAVb3qPkSuY3b2kfU4L+gjIq1xnuI54NYo3cNOLgqAeAARuRtnwT/p1iivDxuBfy25W+d+4KyqHq3pTfV6SkdNtwwWm7mYC/gB75Vcty5Q1QENFnQ9sZmLJsFmLrKAPiLyPVAE/Ieqnm64qOuHzVw8D7wmIv+O8xRGcmM8QBSRd3Gewmtdcr3iRcALQFX/jvP6RSLwM3ABGGVruY0wV4ZhGEYlzC9tDcMwmghT8A3DMJoIU/ANwzCaCFPwDcMwmghT8A3DMJoIU/ANwzCaCFPwDcMwmghT8A3DMJqI/wdNFHuvzFba7QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, Math, Latex\n",
    "import numpy as np\n",
    "import scipy\n",
    "import seaborn\n",
    "import matplotlib.pyplot as plt\n",
    "print(\"suppose we want to sample from a Beta distribution:\")\n",
    "display(Math(r'Beta(x|\\alpha,\\beta)=\\frac{\\Gamma(\\alpha+\\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)}'\n",
    "             r'x^{\\alpha-1}(1-x)^{\\beta-1}'))\n",
    "\n",
    "alpha=2\n",
    "beta=2\n",
    "M=max_pdf=1.5 # The maximum value of the pdf of Beta(x|2,2) is 1.5, so we can use M=1.5 and q(x)=U(0,1)\n",
    "np.random.seed(42)\n",
    "\n",
    "n_samples=1000\n",
    "samples_before_reject=np.random.rand(int(n_samples*M))\n",
    "prob=np.random.rand(int(n_samples*M))\n",
    "rv=scipy.stats.beta(a=2,b=2)\n",
    "retained_samples=[prob<rv.pdf(samples_before_reject)/M]   #This is the reject step\n",
    "\n",
    "final_samples=samples_before_reject[retained_samples]\n",
    "plt.figure(figsize=(6,6))\n",
    "seaborn.distplot(final_samples,\n",
    "                hist=False,\n",
    "                color='r',\n",
    "                label=\"Reject sample\")\n",
    "seaborn.distplot(scipy.stats.beta(a=2,b=2).rvs(n_samples),\n",
    "                hist=False,\n",
    "                color='b',\n",
    "                label='Scipy sample method')\n",
    "plt.xlim([0,1])\n",
    "plt.legend(fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.4  Adaptive rejection sampling\n",
    "We now discuss a method that can automatically come up with a tight upper envelope $q(x)$ to any log concave density $p(x)$. The idea is to upper bound the log density with a piecewise linear function.We choose the initial locations for the pieces based on a fixed grid over the support of the distribution. We then evaluate the gradient of the log density at these locations, and make the lines be tangent at these points.\n",
    "\n",
    "Since the log of the envelope is piecewise linear, the envelope itself is piecewise exponenetial\n",
    "$$\n",
    "q(x)=M_i\\lambda_iexp\\left(-\\lambda_i(x-x_{i-1})\\right),\\quad x \\in[x_{i-1},x_i]\n",
    "$$\n",
    "where x i are the grid points.\n",
    "![1.png](imgs/1.png)\n",
    "\n",
    "## 2.5 Rejection sampling in high dimensions\n",
    "It's clear that we want to make our proposal $q(x)$ as close as possible to the target distribution $p(x)$, while still being an upper bound. But this is quite hard to achieve, especially in high dimensions. To see this, consider sampling from $p(x) = \\mathcal{N} (0, \\sigma_p^2\\mathbf{I})$ using as a proposal $q(x) = \\mathcal{N} (0, \\sigma_q^2\\mathbf{I})$. Obviously we must have $\\sigma_q^2 \\geq \\sigma_p^2$ in order to be an upper bound.In $D$ dimensions, the optimal value is given by $M=\\left(\\frac{\\sigma_q}{\\sigma_p}\\right)^2$. The acceptance rate is $\\frac{1}{M}$, which decreases exponentially fast with dimension.\n",
    "\n",
    "# 3. Importance sampling\n",
    "When we approximate integrals of the form\n",
    "$$\n",
    "\\mathit{I}=\\mathbb{E}\\left[f\\right]=\\int f(x)p(x)dx\n",
    "$$\n",
    "using\n",
    "$$\n",
    "\\hat{\\mathit{I}}=\\frac{1}{S}\\sum_{s=1}^S f(x_s)\n",
    "$$\n",
    "The finite sum approximation to the expectation depends on being able to draw samples from the distribution $p(x)$. Suppose, however, that it is impractical to sample directly from $p(z)$. Furthermore, in many applications  f(x) is nearly zero outside a region $A$ for which $P(x \\in A)$ is small. The set $A$ may have small volume, or it may be in the tail of the distribution. A plain sampling algorithm might lead to high variance because we don't put more attention on $A$.\n",
    "## 3.1 Basic idea\n",
    "It is clear intuitively that we must get some samples from the interesting or important region. We do this by sampling from a distribution that overweights the important region, hence the name **importance sampling**. Having oversampled the important region, we have to adjust our estimate somehow to account for having sampled from this other distribution\n",
    "\n",
    "Importance sampling samples from any proposal $q(x)$. It then uses these samples to estimate the integral as follows\n",
    "\\begin{align}\n",
    "\\mathbb{E}\\left[f\\right]  &=\\int f(x)\\frac{p(x)}{q(x)}q(x)dx  \\\\\n",
    "                          &=\\frac{1}{S}\\sum_{s=1}^S \\frac{p(x^s)}{q(x^s)} f(x^s) \\\\\n",
    "                          &=\\frac{1}{S}\\sum_{s=1}^S w_s\\, f(x^s) \\\\\n",
    "\\end{align}\n",
    "where \n",
    "$$\n",
    "w_s=\\frac{p(x^s)}{q(x^s)}\n",
    "$$\n",
    "are called **importance weights**.Note that, unlike rejection sampling, all of the samples generated are retained but with different weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The importance sampling method result of the variance of a laplace distribution:1.787035\n",
      "The actual variance of a laplace distribution:2.000000\n"
     ]
    }
   ],
   "source": [
    "# Below we will use importance sampling method to calculate the variance of a laplace distribution\n",
    "import scipy\n",
    "import numpy as np\n",
    "\n",
    "n_samples=200\n",
    "np.random.seed(123)\n",
    "samples=scipy.stats.norm.rvs(size=n_samples)\n",
    "samples_variance=np.square(samples)\n",
    "weights=scipy.stats.laplace.pdf(samples)/scipy.stats.norm.pdf(samples)\n",
    "print(\"The importance sampling method result of the variance of a laplace distribution:%4f\"\n",
    "      %(np.mean(samples_variance*weights)))\n",
    "print(\"The actual variance of a laplace distribution:%4f\"%(scipy.stats.laplace.stats(moments='v')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 Handing unnormalized distribution\n",
    "It is frequently the case that we can evaluate the unnormalized target distribution $\\tilde{p}(x)$, but not its normalization constant $Z_p$. We may also want to use an unnormalized proposal $\\tilde{q}(x)$ with possible unknown normalization constant $Z_q$. We can do this as follows.First we evaluate\n",
    "$$\n",
    "\\mathbb{E}\\left[f\\right]=\\frac{Z_q}{Z_p}\\int f(x)\\frac{\\tilde{p}(x)}{\\tilde{q}(x)}q(x)dx \\approx \\frac{Z_q}{Z_p}\\frac{1}{S}\\sum_{s=1}^S\\tilde{w}_sf(x^s)\n",
    "$$\n",
    "where \n",
    "$$\n",
    "\\tilde{w}^s=\\frac{\\tilde{p}(x^s)}{\\tilde{q}(x^s)}\n",
    "$$\n",
    "is the unnormalized importance weight. We could use the same set of samples to evaluate the ratio $\\frac{Z_p}{Z_q}$ as follows.\n",
    "\\begin{align}\n",
    "\\frac{Z_p}{Z_q}  &=\\int \\frac{\\tilde{p}(x)}{Z_q}dx  \\\\\n",
    "                 &=\\int \\frac{\\tilde{p}(x)}{\\tilde{q}(x)}q(x)dx \\\\\n",
    "                 & \\approx \\frac{1}{S}\\sum_{s=1}^S \\tilde{w}^s\n",
    "\\end{align}\n",
    "Hence\n",
    "\\begin{align}\n",
    "\\mathbb{E}\\left[f\\right] &\\approx \\frac{Z_q}{Z_p}\\frac{1}{S}\\sum_{s=1}^S\\tilde{w}_sf(x^s)\\\\\n",
    "                         &\\approx \\frac{\\sum_{s=1}^S\\tilde{w}_sf(x^s)}{\\sum_{s=1}^S \\tilde{w}^s} \\\\\n",
    "                         &=\\sum_{s=1}^S \\bar{w}^s f(x^s)\n",
    "\\end{align}\n",
    "where\n",
    "$$\n",
    "\\bar{w}^s=\\frac{\\tilde{w}_s}{\\sum_{s=1}^S \\tilde{w}^s}\n",
    "$$"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
